PWGLF/STRANGENESS/Cascades/AliAnalysisTaskCheckPerformanceCascadePbPb.cxx

   1 /***************************************************************          *
   2  *  Authors : Antonin Maire, Boris Hippolyte
   3  * Contributors are mentioned in the code where appropriate.              *
   4  *                                                                        *
   5  * Permission to use, copy, modify and distribute this software and its   *
   6  * documentation strictly for non-commercial purposes is hereby granted   *
   7  * without fee, provided that the above copyright notice appears in all   *
   8  * copies and that both the copyright notice and this permission notice   *
   9  * appear in the supporting documentation. The authors make no claims     *
  10  * about the suitability of this software for any purpose. It is          *
  11  * provided "as is" without express or implied warranty.                  *
  12  **************************************************************************/
  13
  14 //-----------------------------------------------------------------
  15 //            AliAnalysisTaskCheckPerformanceCascadePbPb class
  16 //            This task is for a performance study of cascade identification in PbPb.
  17 //            It works with MC info and ESD/AOD.
  18 //            Origin   : AliAnalysisTaskCheckPerformanceCascade class by A. Maire Nov2010, antonin.maire@ires.in2p3.fr
  19 //            Modified for PbPb analysis: M. Nicassio Feb2011, maria.nicassio@ba.infn.it:
  20 //                        - physics selection moved to the run.C macro
  21 //                        - added centrality selection and possibility to select events in nTracks ranges
  22 //                        - added new histograms
  23 //                        - modified binning of some histograms and containers
  24 //                        - flag to enable CF container usage
  25 //                        - check in the destructor for CAF usage
  26 //                        - flag for acceptance cut in the MC part
  27 //                        - in the MC particle selection IsPhysicalPrimary added and number of particles taken as appropriate for HIJING
  28 //                          (however for cascades one gets the same if runs on Nprimaries in the stack and does not check IsPhysicalPrimary)
  29 //                        - automatic settings for PID
  30 //                        - selection of injected cascades and HIJING cascades (kind of "bug" in method IsFromBGEvent())
  31 //                        - added proper time histograms for cascades and lambdas
  32 //                        - cos of PA V0 wrt Xi vertex and not primary vertex
  33 //                        - distance xi-V0 added in the container
  34 //                        - AOD analysis developed (January 2012)
  35 //-----------------------------------------------------------------
  36
  37
  38 #include <Riostream.h>
  39
  40 #include "TList.h"
  41 #include "TFile.h"
  42 #include "TH1F.h"
  43 #include "TH2F.h"
  44 #include "TH3F.h"
  45 #include "TVector3.h"
  46 #include "TCanvas.h"
  47 #include "TParticle.h"
  48 #include "TMath.h"
  49
  50 #include "AliLog.h"
  51 #include "AliHeader.h"
  52 #include "AliMCEvent.h"
  53 #include "AliStack.h"
  54 #include "AliMultiplicity.h"
  55 #include "AliInputEventHandler.h"
  56 #include "AliAnalysisManager.h"
  57
  58 #include "AliCentrality.h"
  59
  60 #include "AliCFContainer.h"
  61
  62 #include "AliESDVZERO.h"
  63
  64 #include "AliGenEventHeader.h"
  65 #include "AliGenCocktailEventHeader.h"
  66 #include "AliGenHijingEventHeader.h"
  67 #include "AliESDtrackCuts.h"
  68 #include "AliPIDResponse.h"
  69 //#include "AliV0vertexer.h"
  70 //#include "AliCascadeVertexer.h"
  71 #include "AliESDEvent.h"
  72 #include "AliESDcascade.h"
  73 #include "AliAODEvent.h"
  74 #include "AliAODMCParticle.h"
  75 #include "AliAnalysisTaskCheckPerformanceCascadePbPb.h"
  76
  77 ClassImp(AliAnalysisTaskCheckPerformanceCascadePbPb)
  78
  79
  80
  81      //_____Dummy constructor________________________________________________________________
  82 AliAnalysisTaskCheckPerformanceCascadePbPb::AliAnalysisTaskCheckPerformanceCascadePbPb()
  83 : AliAnalysisTaskSE(), // <- take care to AliAnalysisTask( empty )
  84   fAnalysisType("ESD"), fESDtrackCuts(0), /*fPaveTextBookKeeping(0),*/
  85
  86     fPIDResponse                   (0),
  87     fkRerunV0CascVertexers         (0),
  88     fkQualityCutZprimVtxPos        (kTRUE),
  89     fkRejectEventPileUp            (kTRUE),
  90     fkQualityCutNoTPConlyPrimVtx   (kTRUE),
  91     fkQualityCutTPCrefit           (kTRUE),
  92     fkQualityCutnTPCcls            (kTRUE),
  93     fMinnTPCcls                    (0),
  94     fkExtraSelections              (0),
  95     fCentrLowLim(0),    fCentrUpLim(0), fCentrEstimator(0),
  96     fVtxRange                      (0),
  97     fApplyAccCut                   (0),
  98
  99
 100         // - Cascade part initialisation
 101     fListHistCascade(0),
 102
 103     // Events in centrality bins
 104     fHistEvtsInCentralityBinsvsNtracks(0),
 105     fHistBestVtxX(0),
 106     fHistBestVtxY(0),
 107     fHistBestVtxZ(0),
 108
 109     // Cascade multiplicity histos
 110
 111     fHistnXiPlusPerEvTot(0),
 112     fHistnXiMinusPerEvTot(0),
 113     fHistnOmegaPlusPerEvTot(0),
 114     fHistnOmegaMinusPerEvTot(0),
 115
 116     fHistnXiPlusPerEv(0),
 117     fHistnXiMinusPerEv(0),
 118     fHistnOmegaPlusPerEv(0),
 119     fHistnOmegaMinusPerEv(0),
 120
 121     fHistnAssoXiMinus(0),
 122     fHistnAssoXiPlus(0),
 123     fHistnAssoOmegaMinus(0),
 124     fHistnAssoOmegaPlus(0),
 125
 126
 127     fHistMCTrackMultiplicity(0),
 128        // - Resolution of the multiplicity estimator
 129     f2dHistRecoMultVsMCMult(0),
 130
 131     fHistEtaGenProton(0),
 132     fHistEtaGenAntiProton(0),
 133
 134    // Xi-
 135    fHistEtaGenCascXiMinus(0),
 136    f3dHistGenPtVsGenYvsCentXiMinusNat(0),
 137    f3dHistGenPtVsGenYvsNtracksXiMinusNat(0),
 138    f3dHistGenPtVsGenYvsCentXiMinusInj(0),
 139    f3dHistGenPtVsGenYvsNtracksXiMinusInj(0),
 140    f3dHistGenPtVsGenctauvsYXiMinusNat(0),
 141    f3dHistGenPtVsGenctauvsYXiMinusInj(0),
 142
 143     fHistThetaGenCascXiMinusNat(0),
 144     fHistThetaGenCascXiMinusInj(0),
 145     f2dHistGenPtVsGenYFdblXiMinus(0),
 146
 147     fHistThetaLambdaXiMinus(0),
 148     fHistThetaBachXiMinus(0),
 149
 150     fHistThetaMesDghterXiMinus(0),
 151     fHistThetaBarDghterXiMinus(0),
 152
 153     fHistPtBachXiMinus(0),
 154     fHistPtMesDghterXiMinus(0),
 155     fHistPtBarDghterXiMinus(0),
 156
 157
 158    // Xi+
 159    fHistEtaGenCascXiPlus(0),
 160    f3dHistGenPtVsGenYvsCentXiPlusNat(0),
 161    f3dHistGenPtVsGenYvsNtracksXiPlusNat(0),
 162    f3dHistGenPtVsGenYvsCentXiPlusInj(0),
 163    f3dHistGenPtVsGenYvsNtracksXiPlusInj(0),
 164    f3dHistGenPtVsGenctauvsYXiPlusNat(0),
 165    f3dHistGenPtVsGenctauvsYXiPlusInj(0),
 166
 167     fHistThetaGenCascXiPlusNat(0),
 168     fHistThetaGenCascXiPlusInj(0),
 169     f2dHistGenPtVsGenYFdblXiPlus(0),
 170
 171     fHistThetaLambdaXiPlus(0),
 172     fHistThetaBachXiPlus(0),
 173
 174     fHistThetaMesDghterXiPlus(0),
 175     fHistThetaBarDghterXiPlus(0),
 176
 177     fHistPtBachXiPlus(0),
 178     fHistPtMesDghterXiPlus(0),
 179     fHistPtBarDghterXiPlus(0),
 180
 181    // Omega-
 182    fHistEtaGenCascOmegaMinus(0),
 183    f3dHistGenPtVsGenYvsCentOmegaMinusNat(0),
 184    f3dHistGenPtVsGenYvsNtracksOmegaMinusNat(0),
 185    f3dHistGenPtVsGenYvsCentOmegaMinusInj(0),
 186    f3dHistGenPtVsGenYvsNtracksOmegaMinusInj(0),
 187    f3dHistGenPtVsGenctauvsYOmegaMinusNat(0),
 188    f3dHistGenPtVsGenctauvsYOmegaMinusInj(0),
 189
 190     fHistThetaGenCascOmegaMinusNat(0),
 191     fHistThetaGenCascOmegaMinusInj(0),
 192     f2dHistGenPtVsGenYFdblOmegaMinus(0),
 193
 194     fHistThetaLambdaOmegaMinus(0),
 195     fHistThetaBachOmegaMinus(0),
 196
 197     fHistThetaMesDghterOmegaMinus(0),
 198     fHistThetaBarDghterOmegaMinus(0),
 199
 200     fHistPtBachOmegaMinus(0),
 201     fHistPtMesDghterOmegaMinus(0),
 202     fHistPtBarDghterOmegaMinus(0),
 203
 204
 205    // Omega+
 206    fHistEtaGenCascOmegaPlus(0),
 207    f3dHistGenPtVsGenYvsCentOmegaPlusNat(0),
 208    f3dHistGenPtVsGenYvsNtracksOmegaPlusNat(0),
 209    f3dHistGenPtVsGenYvsCentOmegaPlusInj(0),
 210    f3dHistGenPtVsGenYvsNtracksOmegaPlusInj(0),
 211    f3dHistGenPtVsGenctauvsYOmegaPlusNat(0),
 212    f3dHistGenPtVsGenctauvsYOmegaPlusInj(0),
 213
 214     fHistThetaGenCascOmegaPlusNat(0),
 215     fHistThetaGenCascOmegaPlusInj(0),
 216     f2dHistGenPtVsGenYFdblOmegaPlus(0),
 217
 218     fHistThetaLambdaOmegaPlus(0),
 219     fHistThetaBachOmegaPlus(0),
 220
 221     fHistThetaMesDghterOmegaPlus(0),
 222     fHistThetaBarDghterOmegaPlus(0),
 223
 224     fHistPtBachOmegaPlus(0),
 225     fHistPtMesDghterOmegaPlus(0),
 226     fHistPtBarDghterOmegaPlus(0),
 227
 228 // Part 2 - Association to MC
 229
 230     fHistMassXiMinus(0),
 231     fHistMassXiPlus(0),
 232     fHistMassOmegaMinus(0),
 233     fHistMassOmegaPlus(0),
 234
 235         // - Effective mass histos with combined PID
 236     fHistMassWithCombPIDXiMinus(0), fHistMassWithCombPIDXiPlus(0),
 237     fHistMassWithCombPIDOmegaMinus(0), fHistMassWithCombPIDOmegaPlus(0),
 238
 239         // - PID Probability versus MC Pt(bachelor track)
 240     f2dHistPIDprobaKaonVsMCPtBach(0), f2dHistPIDprobaPionVsMCPtBach(0),
 241
 242         // - Effective mass histos with perfect MC PID on the bachelor
 243     fHistMassWithMcPIDXiMinus(0), fHistMassWithMcPIDXiPlus(0),
 244     fHistMassWithMcPIDOmegaMinus(0), fHistMassWithMcPIDOmegaPlus(0),
 245
 246
 247         // - Effective mass histos for the cascade candidates associated with MC
 248     fHistAsMCMassXiMinus(0),
 249     fHistAsMCMassXiPlus(0),
 250     fHistAsMCMassOmegaMinus(0),
 251     fHistAsMCMassOmegaPlus(0),
 252
 253         // - Generated Pt Vs generated y, for the cascade candidates associated with MC + Info Comb. PID
 254     f2dHistAsMCandCombPIDGenPtVsGenYXiMinus(0),
 255     f2dHistAsMCandCombPIDGenPtVsGenYXiPlus(0),
 256     f2dHistAsMCandCombPIDGenPtVsGenYOmegaMinus(0),
 257     f2dHistAsMCandCombPIDGenPtVsGenYOmegaPlus(0),
 258
 259         // - Generated Pt Vs generated y, for the cascade candidates associated with MC
 260     f2dHistAsMCGenPtVsGenYXiMinus(0),
 261     f2dHistAsMCGenPtVsGenYXiPlus(0),
 262     f2dHistAsMCGenPtVsGenYOmegaMinus(0),
 263     f2dHistAsMCGenPtVsGenYOmegaPlus(0),
 264
 265         // - Generated Eta of the the cascade candidates associated with MC
 266     fHistAsMCGenEtaXiMinus(0),
 267     fHistAsMCGenEtaXiPlus(0),
 268     fHistAsMCGenEtaOmegaMinus(0),
 269     fHistAsMCGenEtaOmegaPlus(0),
 270
 271         // - Resolution in Pt as function of generated Pt
 272     f2dHistAsMCResPtXiMinus(0),
 273     f2dHistAsMCResPtXiPlus(0),
 274     f2dHistAsMCResPtOmegaMinus(0),
 275     f2dHistAsMCResPtOmegaPlus(0),
 276
 277         // - Resolution in R(2D) as function of generated R
 278     f2dHistAsMCResRXiMinus(0),
 279     f2dHistAsMCResRXiPlus(0),
 280     f2dHistAsMCResROmegaMinus(0),
 281     f2dHistAsMCResROmegaPlus(0),
 282
 283         // - Resolution in phi as function of generated Pt
 284     f2dHistAsMCResPhiXiMinus(0),
 285     f2dHistAsMCResPhiXiPlus(0),
 286     f2dHistAsMCResPhiOmegaMinus(0),
 287     f2dHistAsMCResPhiOmegaPlus(0),
 288
 289        //  - Correlation between proton (antiproton) daughter MC pt and Xi/Omega MC pt (to apply Geat/Fluka correction)
 290     f2dHistAsMCptProtonMCptXiMinus(0),
 291     f2dHistAsMCptAntiprotonMCptXiPlus(0),
 292     f2dHistAsMCptProtonMCptOmegaMinus(0),
 293     f2dHistAsMCptAntiprotonMCptOmegaPlus(0),
 294
 295     fHistV0toXiCosineOfPointingAngle(0),
 296     fHistV0CosineOfPointingAnglevsPtXi(0),
 297     fHistV0CosineOfPointingAnglevsPtOmega(0),
 298
 299     fCFContCascadePIDAsXiMinus(0),
 300     fCFContCascadePIDAsXiPlus(0),
 301     fCFContCascadePIDAsOmegaMinus(0),
 302     fCFContCascadePIDAsOmegaPlus(0),
 303
 304     fCFContAsCascadeCuts(0),
 305
 306     fV0Ampl             (0)
 307
 308 {
 309 // Dummy constructor
 310         for(Int_t iV0selIdx   = 0; iV0selIdx   < 7; iV0selIdx++   ) { fV0Sels          [iV0selIdx   ] = -1.; }
 311         for(Int_t iCascSelIdx = 0; iCascSelIdx < 8; iCascSelIdx++ ) { fCascSels        [iCascSelIdx ] = -1.; }
 312 }
 313
 314
 315
 316
 317 //_____Non-default Constructor________________________________________________________________
 318 AliAnalysisTaskCheckPerformanceCascadePbPb::AliAnalysisTaskCheckPerformanceCascadePbPb(const char *name)
 319   : AliAnalysisTaskSE(name),
 320     fAnalysisType("ESD"), fESDtrackCuts(0), /*fPaveTextBookKeeping(0),*/
 321
 322     fPIDResponse                   (0),
 323     fkRerunV0CascVertexers         (0),
 324     fkQualityCutZprimVtxPos        (kTRUE),
 325     fkRejectEventPileUp            (kTRUE),
 326     fkQualityCutNoTPConlyPrimVtx   (kTRUE),
 327     fkQualityCutTPCrefit           (kTRUE),
 328     fkQualityCutnTPCcls            (kTRUE),
 329     fMinnTPCcls                    (0),
 330     fkExtraSelections              (0),
 331     fCentrLowLim(0), fCentrUpLim(0), fCentrEstimator(0),
 332     fVtxRange                      (0),
 333     fApplyAccCut                   (0),
 334
 335
 336         // - Cascade part initialisation
 337     fListHistCascade(0),
 338
 339     // Events in centraity bins
 340     fHistEvtsInCentralityBinsvsNtracks(0),
 341     fHistBestVtxX(0),
 342     fHistBestVtxY(0),
 343     fHistBestVtxZ(0),
 344
 345     // Cascade multiplicity histos
 346
 347     fHistnXiPlusPerEvTot(0),
 348     fHistnXiMinusPerEvTot(0),
 349     fHistnOmegaPlusPerEvTot(0),
 350     fHistnOmegaMinusPerEvTot(0),
 351
 352     fHistnXiPlusPerEv(0),
 353     fHistnXiMinusPerEv(0),
 354     fHistnOmegaPlusPerEv(0),
 355     fHistnOmegaMinusPerEv(0),
 356
 357     fHistnAssoXiMinus(0),
 358     fHistnAssoXiPlus(0),
 359     fHistnAssoOmegaMinus(0),
 360     fHistnAssoOmegaPlus(0),
 361
 362     fHistMCTrackMultiplicity(0),
 363        // - Resolution of the multiplicity estimator
 364     f2dHistRecoMultVsMCMult(0),
 365
 366     fHistEtaGenProton(0),
 367     fHistEtaGenAntiProton(0),
 368
 369 // Xi-
 370    fHistEtaGenCascXiMinus(0),
 371    f3dHistGenPtVsGenYvsCentXiMinusNat(0),
 372    f3dHistGenPtVsGenYvsNtracksXiMinusNat(0),
 373    f3dHistGenPtVsGenYvsCentXiMinusInj(0),
 374    f3dHistGenPtVsGenYvsNtracksXiMinusInj(0),
 375    f3dHistGenPtVsGenctauvsYXiMinusNat(0),
 376    f3dHistGenPtVsGenctauvsYXiMinusInj(0),
 377
 378     fHistThetaGenCascXiMinusNat(0),
 379     fHistThetaGenCascXiMinusInj(0),
 380     f2dHistGenPtVsGenYFdblXiMinus(0),
 381
 382     fHistThetaLambdaXiMinus(0),
 383     fHistThetaBachXiMinus(0),
 384
 385     fHistThetaMesDghterXiMinus(0),
 386     fHistThetaBarDghterXiMinus(0),
 387
 388     fHistPtBachXiMinus(0),
 389     fHistPtMesDghterXiMinus(0),
 390     fHistPtBarDghterXiMinus(0),
 391
 392
 393    // Xi+
 394    fHistEtaGenCascXiPlus(0),
 395   f3dHistGenPtVsGenYvsCentXiPlusNat(0),
 396   f3dHistGenPtVsGenYvsNtracksXiPlusNat(0),
 397   f3dHistGenPtVsGenYvsCentXiPlusInj(0),
 398   f3dHistGenPtVsGenYvsNtracksXiPlusInj(0),
 399   f3dHistGenPtVsGenctauvsYXiPlusNat(0),
 400   f3dHistGenPtVsGenctauvsYXiPlusInj(0),
 401
 402     fHistThetaGenCascXiPlusNat(0),
 403     fHistThetaGenCascXiPlusInj(0),
 404     f2dHistGenPtVsGenYFdblXiPlus(0),
 405
 406     fHistThetaLambdaXiPlus(0),
 407     fHistThetaBachXiPlus(0),
 408
 409     fHistThetaMesDghterXiPlus(0),
 410     fHistThetaBarDghterXiPlus(0),
 411
 412     fHistPtBachXiPlus(0),
 413     fHistPtMesDghterXiPlus(0),
 414     fHistPtBarDghterXiPlus(0),
 415
 416    // Omega-
 417    fHistEtaGenCascOmegaMinus(0),
 418    f3dHistGenPtVsGenYvsCentOmegaMinusNat(0),
 419    f3dHistGenPtVsGenYvsNtracksOmegaMinusNat(0),
 420    f3dHistGenPtVsGenYvsCentOmegaMinusInj(0),
 421    f3dHistGenPtVsGenYvsNtracksOmegaMinusInj(0),
 422    f3dHistGenPtVsGenctauvsYOmegaMinusNat(0),
 423    f3dHistGenPtVsGenctauvsYOmegaMinusInj(0),
 424
 425     fHistThetaGenCascOmegaMinusNat(0),
 426     fHistThetaGenCascOmegaMinusInj(0),
 427     f2dHistGenPtVsGenYFdblOmegaMinus(0),
 428
 429     fHistThetaLambdaOmegaMinus(0),
 430     fHistThetaBachOmegaMinus(0),
 431
 432     fHistThetaMesDghterOmegaMinus(0),
 433     fHistThetaBarDghterOmegaMinus(0),
 434
 435     fHistPtBachOmegaMinus(0),
 436     fHistPtMesDghterOmegaMinus(0),
 437     fHistPtBarDghterOmegaMinus(0),
 438
 439
 440    // Omega+
 441    fHistEtaGenCascOmegaPlus(0),
 442    f3dHistGenPtVsGenYvsCentOmegaPlusNat(0),
 443    f3dHistGenPtVsGenYvsNtracksOmegaPlusNat(0),
 444    f3dHistGenPtVsGenYvsCentOmegaPlusInj(0),
 445    f3dHistGenPtVsGenYvsNtracksOmegaPlusInj(0),
 446    f3dHistGenPtVsGenctauvsYOmegaPlusNat(0),
 447    f3dHistGenPtVsGenctauvsYOmegaPlusInj(0),
 448
 449     fHistThetaGenCascOmegaPlusNat(0),
 450     fHistThetaGenCascOmegaPlusInj(0),
 451     f2dHistGenPtVsGenYFdblOmegaPlus(0),
 452
 453     fHistThetaLambdaOmegaPlus(0),
 454     fHistThetaBachOmegaPlus(0),
 455
 456     fHistThetaMesDghterOmegaPlus(0),
 457     fHistThetaBarDghterOmegaPlus(0),
 458
 459     fHistPtBachOmegaPlus(0),
 460     fHistPtMesDghterOmegaPlus(0),
 461     fHistPtBarDghterOmegaPlus(0),
 462
 463 // Part 2 - Association to MC
 464
 465     fHistMassXiMinus(0),
 466     fHistMassXiPlus(0),
 467     fHistMassOmegaMinus(0),
 468     fHistMassOmegaPlus(0),
 469
 470         // - Effective mass histos with combined PID
 471     fHistMassWithCombPIDXiMinus(0), fHistMassWithCombPIDXiPlus(0),
 472     fHistMassWithCombPIDOmegaMinus(0), fHistMassWithCombPIDOmegaPlus(0),
 473
 474         // - PID Probability versus MC Pt(bachelor track)
 475     f2dHistPIDprobaKaonVsMCPtBach(0), f2dHistPIDprobaPionVsMCPtBach(0),
 476
 477         // - Effective mass histos with perfect MC PID on the bachelor
 478     fHistMassWithMcPIDXiMinus(0), fHistMassWithMcPIDXiPlus(0),
 479     fHistMassWithMcPIDOmegaMinus(0), fHistMassWithMcPIDOmegaPlus(0),
 480
 481
 482         // - Effective mass histos for the cascade candidates associated with MC
 483     fHistAsMCMassXiMinus(0),
 484     fHistAsMCMassXiPlus(0),
 485     fHistAsMCMassOmegaMinus(0),
 486     fHistAsMCMassOmegaPlus(0),
 487
 488         // - Generated Pt Vs generated y, for the cascade candidates associated with MC + Info Comb. PID
 489     f2dHistAsMCandCombPIDGenPtVsGenYXiMinus(0),
 490     f2dHistAsMCandCombPIDGenPtVsGenYXiPlus(0),
 491     f2dHistAsMCandCombPIDGenPtVsGenYOmegaMinus(0),
 492     f2dHistAsMCandCombPIDGenPtVsGenYOmegaPlus(0),
 493
 494         // - Generated Pt Vs generated y, for the cascade candidates associated with MC
 495     f2dHistAsMCGenPtVsGenYXiMinus(0),
 496     f2dHistAsMCGenPtVsGenYXiPlus(0),
 497     f2dHistAsMCGenPtVsGenYOmegaMinus(0),
 498     f2dHistAsMCGenPtVsGenYOmegaPlus(0),
 499
 500         // - Generated Eta of the the cascade candidates associated with MC
 501     fHistAsMCGenEtaXiMinus(0),
 502     fHistAsMCGenEtaXiPlus(0),
 503     fHistAsMCGenEtaOmegaMinus(0),
 504     fHistAsMCGenEtaOmegaPlus(0),
 505
 506         // - Resolution in Pt as function of generated Pt
 507     f2dHistAsMCResPtXiMinus(0),
 508     f2dHistAsMCResPtXiPlus(0),
 509     f2dHistAsMCResPtOmegaMinus(0),
 510     f2dHistAsMCResPtOmegaPlus(0),
 511
 512         // - Resolution in R(2D) as function of generated R
 513     f2dHistAsMCResRXiMinus(0),
 514     f2dHistAsMCResRXiPlus(0),
 515     f2dHistAsMCResROmegaMinus(0),
 516     f2dHistAsMCResROmegaPlus(0),
 517
 518        // - Resolution in phi as function of generated Pt
 519     f2dHistAsMCResPhiXiMinus(0),
 520     f2dHistAsMCResPhiXiPlus(0),
 521     f2dHistAsMCResPhiOmegaMinus(0),
 522     f2dHistAsMCResPhiOmegaPlus(0),
 523
 524        //  - Correlation between proton (antiproton) daughter MC pt and Xi/Omega MC pt (to apply Geat/Fluka correction)
 525     f2dHistAsMCptProtonMCptXiMinus(0),
 526     f2dHistAsMCptAntiprotonMCptXiPlus(0),
 527     f2dHistAsMCptProtonMCptOmegaMinus(0),
 528     f2dHistAsMCptAntiprotonMCptOmegaPlus(0),
 529
 530     fHistV0toXiCosineOfPointingAngle(0),
 531     fHistV0CosineOfPointingAnglevsPtXi(0),
 532     fHistV0CosineOfPointingAnglevsPtOmega(0),
 533
 534     fCFContCascadePIDAsXiMinus(0),
 535     fCFContCascadePIDAsXiPlus(0),
 536     fCFContCascadePIDAsOmegaMinus(0),
 537     fCFContCascadePIDAsOmegaPlus(0),
 538
 539     fCFContAsCascadeCuts(0),
 540
 541     fV0Ampl(0)
 542
 543 {
 544   // Constructor
 545
 546   // Define input and output slots here
 547   // Input slot #0 works with a TChain
 548   // Output slot #1 writes into a TList container (cascade)
 549
 550
 551         // PbPb default cuts
 552
 553         fV0Sels[0] =  33.  ;  // max allowed chi2
 554         fV0Sels[1] =   0.1; // min allowed impact parameter for the 1st daughter
 555         fV0Sels[2] =   0.1; // min allowed impact parameter for the 2nd daughter
 556         fV0Sels[3] =   1.0 ;  // max allowed DCA between the daughter tracks
 557         fV0Sels[4] =   0.998 ;  // min allowed cosine of V0's pointing angle
 558         fV0Sels[5] =   0.9;  // min radius of the fiducial volume
 559         fV0Sels[6] = 100.  ;  // max radius of the fiducial volume
 560
 561         fCascSels[0] =  33.   ;  // max allowed chi2
 562         fCascSels[1] =   0.05;  // min allowed V0 impact parameter
 563         fCascSels[2] =   0.008;  // "window" around the Lambda mass
 564         fCascSels[3] =   0.03;  // min allowed bachelor's impact parameter
 565         fCascSels[4] =   0.3  ;  // max allowed DCA between the V0 and the bachelor
 566         fCascSels[5] =   0.999;  // min allowed cosine of the cascade pointing angle
 567         fCascSels[6] =   0.9  ;  // min radius of the fiducial volume
 568         fCascSels[7] = 100.   ;  // max radius of the fiducial volume
 569
 570
 571   DefineOutput(1, TList::Class());
 572
 573 }
 574
 575
 576 AliAnalysisTaskCheckPerformanceCascadePbPb::~AliAnalysisTaskCheckPerformanceCascadePbPb()
 577 {
 578   //
 579   // Destructor
 580   //
 581
 582   // For all TH1, 2, 3 HnSparse and CFContainer are in the fListCascade TList.
 583   // They will be deleted when fListCascade is deleted by the TSelector dtor
 584   // Because of TList::SetOwner()
 585
 586   if (fListHistCascade && !AliAnalysisManager::GetAnalysisManager()->IsProofMode())      { delete fListHistCascade;     fListHistCascade = 0x0;  }
 587   if (fESDtrackCuts)         { delete fESDtrackCuts;        fESDtrackCuts = 0x0; }
 588   /*if (fPaveTextBookKeeping)  { delete fPaveTextBookKeeping; fPaveTextBookKeeping = 0x0; } // fPaveTextBookKeeping is not stored into the TList*/
 589 }
 590
 591
 592 //________________________________________________________________________
 593 void AliAnalysisTaskCheckPerformanceCascadePbPb::UserCreateOutputObjects()
 594 {
 595   // Create histograms
 596   // Called once
 597
 598    // Option for AliLog
 599         AliLog::SetGlobalLogLevel(AliLog::kError);
 600         // to suppress the extensive info prompted by a run with MC
 601
 602    // Definition of the output datamembers
 603    fListHistCascade = new TList();
 604    fListHistCascade->SetOwner(); // See http://root.cern.ch/root/html/TCollection.html#TCollection:SetOwner
 605
 606    // New PID object
 607    AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 608    AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 609    fPIDResponse = inputHandler->GetPIDResponse();
 610
 611    // Only used to get the number of primary reconstructed tracks
 612    if (! fESDtrackCuts ){
 613      fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE); // Std definition of primary (see kTRUE argument) tracks for 2010
 614 //     fESDtrackCuts->SetEtaRange(-0.8,+0.8);
 615 //     fESDtrackCuts->SetPtRange(0.15, 1e10);
 616      Printf("CheckCascade - ESDtrackCuts set up to 2010 std ITS-TPC cuts...");
 617    }
 618
 619
 620 /*
 621 if( !fPaveTextBookKeeping){
 622         fPaveTextBookKeeping = new TPaveText(0.1, 0.1, 0.9, 0.9,"NDC");
 623         fPaveTextBookKeeping->SetName("fPaveTextBookKeeping");
 624         fPaveTextBookKeeping->SetBorderSize(0);
 625         fPaveTextBookKeeping->SetTextAlign(12);
 626         fPaveTextBookKeeping->SetFillColor(kWhite);
 627         fPaveTextBookKeeping->SetTextFont(42);        // regular Arial or Helvetica,
 628         fPaveTextBookKeeping->SetTextColor(kGray+3);
 629
 630
 631         fPaveTextBookKeeping->AddText( "Task CHECK PERFORMANCE CASCADE analysis" );
 632         fPaveTextBookKeeping->AddText("- - - - - - - - - - - ");
 633         fPaveTextBookKeeping->AddText( Form("AnalysisType : %s ", fAnalysisType.Data() ));
 634         fPaveTextBookKeeping->AddText("Colliding system : A-A collisions ");
 635
 636         fPaveTextBookKeeping->AddText("- - - - - - - - - - - ");
 637
 638         if(fkRerunV0CascVertexers){
 639                 fPaveTextBookKeeping->AddText("A.1. With V0 vertexer : ");
 640                 fPaveTextBookKeeping->AddText( Form("  - V0 #chi^{2} _________________ <  %.3f ",               fV0Sels[0] ));
 641                 fPaveTextBookKeeping->AddText( Form("  - DCA(prim. Vtx/ 1^{st} daughter) ___ >  %.3f     cm ",  fV0Sels[1] ));
 642                 fPaveTextBookKeeping->AddText( Form("  - DCA(prim. Vtx/ 2^{nd} daughter) __  >  %.3f     cm",   fV0Sels[2] ));
 643                 fPaveTextBookKeeping->AddText( Form("  - DCA between V0 daughters ___ <  %.3f      cm",         fV0Sels[3] ));
 644                 fPaveTextBookKeeping->AddText( Form("  - cos(V0 pointing angle) _______ >  %.3f ",              fV0Sels[4] ));
 645                 fPaveTextBookKeeping->AddText( Form("  - R_{transv}(V0 decay) ________ >  %.3f             cm", fV0Sels[5] ));
 646                 fPaveTextBookKeeping->AddText( Form("  - R_{transv}(V0 decay) ________ <  %.3f         cm",     fV0Sels[6] ));
 647
 648                 fPaveTextBookKeeping->AddText(" ");
 649
 650                 fPaveTextBookKeeping->AddText("A.2. With Casc. vertexer : ");
 651                 fPaveTextBookKeeping->AddText( Form("  - Casc. #chi^{2} ______________  <  %.3f ",                               fCascSels[0] ));
 652                 fPaveTextBookKeeping->AddText( Form("  - DCA(prim. Vtx/ V0) _________ >  %.3f    cm",                            fCascSels[1] ));
 653                 fPaveTextBookKeeping->AddText( Form("  - | M_{#Lambda}(reco) - M_{#Lambda}(pdg) | _______ <  %.3f    GeV/c^{2}", fCascSels[2] ));
 654                 fPaveTextBookKeeping->AddText( Form("  - DCA(prim. Vtx/ Bach) _______ >  %.3f    cm",                            fCascSels[3] ));
 655                 fPaveTextBookKeeping->AddText( Form("  - DCA between Bach/ #Lambda ______ <  %.3f    cm",                        fCascSels[4] ));
 656                 fPaveTextBookKeeping->AddText( Form("  - cos(Casc. pointing angle) ____ >  %.3f ",                               fCascSels[5] ));
 657                 fPaveTextBookKeeping->AddText( Form("  - R_{transv}(Casc. decay) ______ >  %.3f         cm",                     fCascSels[6] ));
 658                 fPaveTextBookKeeping->AddText( Form("  - R_{transv}(Casc. decay) ______ <  %.3f     cm",                         fCascSels[7] ));
 659         }
 660         else{   fPaveTextBookKeeping->AddText("A. No rerunning of the V0/Casc. vertexers ... See std cuts in (AliRoot+Rec.C) used for this prod. cycle");}
 661
 662         fPaveTextBookKeeping->AddText("- - - - - - - - - - - ");
 663
 664         if(fkQualityCutZprimVtxPos)      fPaveTextBookKeeping->AddText("B. Quality Cut(prim. Vtx z-Pos)    = ON  ");
 665         else                             fPaveTextBookKeeping->AddText("B. Quality Cut(prim. Vtx z-Pos)    = Off ");
 666         if(fkQualityCutNoTPConlyPrimVtx) fPaveTextBookKeeping->AddText("C. Quality Cut(No TPC prim. vtx) = ON  ");
 667         else                             fPaveTextBookKeeping->AddText("C. Quality Cut(No TPC prim. vtx) = Off ");
 668         if(fkQualityCutTPCrefit)         fPaveTextBookKeeping->AddText("D. Quality Cut(TPCrefit)               = ON  ");
 669         else                             fPaveTextBookKeeping->AddText("D. Quality Cut(TPCrefit)               = Off ");
 670         if(fkQualityCutnTPCcls)          fPaveTextBookKeeping->AddText("E. Quality Cut(n TPC clusters)   = ON  ");
 671         else                             fPaveTextBookKeeping->AddText("E. Quality Cut(n TPC clusters)   = Off ");
 672         if(fkExtraSelections)            fPaveTextBookKeeping->AddText("F. Extra Analysis Selections         = ON  ");
 673         else                             fPaveTextBookKeeping->AddText("F. Extra Analysis Selections         = Off ");
 674
 675         fPaveTextBookKeeping->AddText("- - - - - - - - - - - ");
 676
 677         fListHistCascade->Add(fPaveTextBookKeeping);
 678 }
 679 */
 680
 681   // - General
 682   Double_t ptBinLimits[101];
 683   for (Int_t iptbin = 0; iptbin<101; ++iptbin) ptBinLimits[iptbin]=iptbin*0.1;
 684   Double_t yBinLimits[111];
 685   for (Int_t iybin = 0; iybin<111; ++iybin) yBinLimits[iybin]=-1.1+iybin*0.02;
 686
 687   // Events in centrality bins
 688   Double_t centBinLimits[12] = {0.,5.,10.,20.,30.,40.,50.,60.,70.,80.,90.,100.};
 689
 690   Double_t ctauBinLimits[112];
 691   for (Int_t ict = 0; ict<112; ++ict) ctauBinLimits[ict] = (Double_t) (ict-1.);
 692
 693   fHistEvtsInCentralityBinsvsNtracks = new TH2F("fHistEvtsInCentralityBinsvsNtracks","",11,centBinLimits,100,0.,4000.);
 694   fListHistCascade->Add(fHistEvtsInCentralityBinsvsNtracks);
 695
 696   fHistBestVtxX = new TH1F("fHistBestVtxX","",100,-0.1,0.1);
 697   fListHistCascade->Add(fHistBestVtxX);
 698   fHistBestVtxY = new TH1F("fHistBestVtxY","",100,-1,1);
 699   fListHistCascade->Add(fHistBestVtxY);
 700   fHistBestVtxZ = new TH1F("fHistBestVtxZ","",100,-20.,20.);
 701   fListHistCascade->Add(fHistBestVtxZ);
 702
 703   // Cascade multiplicity distributions
 704
 705   fHistnXiPlusPerEvTot= new TH1F("fHistnXiPlusPerEvTot", "", 100, 0, 100);
 706   fListHistCascade->Add(fHistnXiPlusPerEvTot);
 707   fHistnXiMinusPerEvTot= new TH1F("fHistnXiMinusPerEvTot", "", 100, 0, 100);
 708   fListHistCascade->Add(fHistnXiMinusPerEvTot);
 709   fHistnOmegaPlusPerEvTot = new TH1F("fHistnOmegaPlusPerEvTot", "", 50, 0, 50);
 710   fListHistCascade->Add(fHistnOmegaPlusPerEvTot);
 711   fHistnOmegaMinusPerEvTot= new TH1F("fHistnOmegaMinusPerEvTot", "", 50, 0, 50);
 712   fListHistCascade->Add(fHistnOmegaMinusPerEvTot);
 713
 714   fHistnXiPlusPerEv= new TH1F("fHistnXiPlusPerEv", "", 100, 0, 100);
 715   fListHistCascade->Add(fHistnXiPlusPerEv);
 716   fHistnXiMinusPerEv= new TH1F("fHistnXiMinusPerEv", "", 100, 0, 100);
 717   fListHistCascade->Add(fHistnXiMinusPerEv);
 718   fHistnOmegaPlusPerEv= new TH1F("fHistnOmegaPlusPerEv", "", 50, 0, 50);
 719   fListHistCascade->Add(fHistnOmegaPlusPerEv);
 720   fHistnOmegaMinusPerEv= new TH1F("fHistnOmegaMinusPerEv", "", 50, 0, 50);
 721   fListHistCascade->Add(fHistnOmegaMinusPerEv);
 722
 723   fHistnAssoXiMinus= new TH1F("fHistnAssoXiMinus", "", 100, 0, 100);
 724   fListHistCascade->Add(fHistnAssoXiMinus);
 725   fHistnAssoXiPlus= new TH1F("fHistnAssoXiPlus", "", 100, 0, 100);
 726   fListHistCascade->Add(fHistnAssoXiPlus);
 727   fHistnAssoOmegaMinus= new TH1F("fHistnAssoOmegaMinus", "", 50, 0, 50);
 728   fListHistCascade->Add(fHistnAssoOmegaMinus);
 729   fHistnAssoOmegaPlus= new TH1F("fHistnAssoOmegaPlus", "", 50, 0, 50);
 730   fListHistCascade->Add(fHistnAssoOmegaPlus);
 731
 732   if (!fHistMCTrackMultiplicity) {
 733     fHistMCTrackMultiplicity = new TH1F("fHistMCTrackMultiplicity", "Multiplicity distribution;Number of tracks;Events", 200, 0, 200000.);
 734     fListHistCascade->Add(fHistMCTrackMultiplicity);
 735   }
 736
 737     // - Resolution of the multiplicity estimator
 738   if (! f2dHistRecoMultVsMCMult){
 739     f2dHistRecoMultVsMCMult = new TH2F("f2dHistRecoMultVsMCMult", "Resolution of the multiplicity estimator (prim. tracks); Reco Multiplicity (prim. tracks); MC multiplicity (gen. part.)", 200, 0., 150000., 200, 0., 200000.);
 740     fListHistCascade->Add(f2dHistRecoMultVsMCMult);
 741   }
 742
 743   if (!fHistEtaGenProton) {
 744     fHistEtaGenProton = new TH1F("fHistEtaGenProton", "#eta of any gen. p^{+};#eta;Number of prim. protons", 200, -10, 10);
 745     fListHistCascade->Add(fHistEtaGenProton);
 746   }
 747
 748   if (!fHistEtaGenAntiProton) {
 749      fHistEtaGenAntiProton = new TH1F("fHistEtaGenAntiProton", "#eta of any gen. #bar{p}^{-};#eta;Number of prim. #bar{p}", 200, -10, 10);
 750      fListHistCascade->Add(fHistEtaGenAntiProton);
 751   }
 752
 753
 754   //--------
 755   // I - Xi-
 756   // - Pseudo-Rapidity distribution
 757   if (!fHistEtaGenCascXiMinus) {
 758      fHistEtaGenCascXiMinus = new TH1F("fHistEtaGenCascXiMinus", "#eta of any gen. #Xi^{-};#eta;Number of Casc", 200, -10, 10);
 759      fListHistCascade->Add(fHistEtaGenCascXiMinus);
 760    }
 761
 762   if (!f3dHistGenPtVsGenYvsCentXiMinusNat) {
 763      f3dHistGenPtVsGenYvsCentXiMinusNat = new TH3D("f3dHistGenPtVsGenYvsCentXiMinusNat", "MC P_{t} Vs MC Y of Gen #Xi^{-} ;Pt_{MC} (GeV/c); Y_{MC}", 100, ptBinLimits, 110, yBinLimits, 11, centBinLimits);
 764      fListHistCascade->Add(f3dHistGenPtVsGenYvsCentXiMinusNat);
 765    }
 766   if (!f3dHistGenPtVsGenYvsNtracksXiMinusNat) {
 767     f3dHistGenPtVsGenYvsNtracksXiMinusNat = new TH3D("f3dHistGenPtVsGenYvsNtracksXiMinusNat", "MC P_{t} Vs MC Y of Gen #Xi^{-} ;Pt_{MC} (GeV/c); Y_{MC}", 100, 0., 10., 110, -1.1, 1.1, 100, 0., 4000.);
 768     fListHistCascade->Add(f3dHistGenPtVsGenYvsNtracksXiMinusNat);
 769   }
 770   if (!f3dHistGenPtVsGenYvsCentXiMinusInj) {
 771      f3dHistGenPtVsGenYvsCentXiMinusInj = new TH3D("f3dHistGenPtVsGenYvsCentXiMinusInj", "MC P_{t} Vs MC Y of Gen #Xi^{-} ;Pt_{MC} (GeV/c); Y_{MC}", 100, ptBinLimits, 110, yBinLimits, 11, centBinLimits);
 772      fListHistCascade->Add(f3dHistGenPtVsGenYvsCentXiMinusInj);
 773   }
 774   if (!f3dHistGenPtVsGenYvsNtracksXiMinusInj) {
 775     f3dHistGenPtVsGenYvsNtracksXiMinusInj = new TH3D("f3dHistGenPtVsGenYvsNtracksXiMinusInj", "MC P_{t} Vs MC Y of Gen #Xi^{-} ;Pt_{MC} (GeV/c); Y_{MC}", 100, 0., 10., 110, -1.1, 1.1, 100, 0., 4000.);
 776     fListHistCascade->Add(f3dHistGenPtVsGenYvsNtracksXiMinusInj);
 777   }
 778   if (!f3dHistGenPtVsGenctauvsYXiMinusNat) {
 779      f3dHistGenPtVsGenctauvsYXiMinusNat = new TH3D("f3dHistGenPtVsGenctauvsYXiMinusNat", "MC P_{t} Vs MC ctau Vs Y of Gen #Xi^{-} ", 100, ptBinLimits, 111, ctauBinLimits, 110, yBinLimits);
 780      fListHistCascade->Add(f3dHistGenPtVsGenctauvsYXiMinusNat);
 781   }
 782
 783   if (!f3dHistGenPtVsGenctauvsYXiMinusInj) {
 784      f3dHistGenPtVsGenctauvsYXiMinusInj = new TH3D("f3dHistGenPtVsGenctauvsYXiMinusInj", "MC P_{t} Vs MC ctau Vs Y of Gen #Xi^{-} ", 100, ptBinLimits, 111, ctauBinLimits, 110, yBinLimits);
 785      fListHistCascade->Add(f3dHistGenPtVsGenctauvsYXiMinusInj);
 786   }
 787
 788
 789   // - Info at the generation level of multi-strange particle
 790
 791   if (!fHistThetaGenCascXiMinusNat) {
 792      fHistThetaGenCascXiMinusNat = new TH1F("fHistThetaGenCascXiMinusNat", "#theta of gen. #Xi^{-};#theta;Number of Casc.", 200, -10, 190);
 793      fListHistCascade->Add(fHistThetaGenCascXiMinusNat);
 794   }
 795
 796   if (!fHistThetaGenCascXiMinusInj) {
 797      fHistThetaGenCascXiMinusInj = new TH1F("fHistThetaGenCascXiMinusInj", "#theta of injected. #Xi^{-};#theta;Number of Casc.", 200, -10, 190);
 798      fListHistCascade->Add(fHistThetaGenCascXiMinusInj);
 799   }
 800
 801
 802   if (!f2dHistGenPtVsGenYFdblXiMinus) {
 803      f2dHistGenPtVsGenYFdblXiMinus = new TH2D("f2dHistGenPtVsGenYFdblXiMinus", "MC P_{t} Vs MC Y of findable Gen #Xi^{-}; Pt_{MC} (GeV/c); Y_{MC}", 100, 0., 10., 220, -1.1, 1.1);
 804      fListHistCascade->Add(f2dHistGenPtVsGenYFdblXiMinus);
 805   }
 806
 807                 // - Theta distribution the daughters (control plots)
 808
 809   if (!fHistThetaLambdaXiMinus) {
 810      fHistThetaLambdaXiMinus = new TH1F("fHistThetaLambdaXiMinus", "#theta of gen. #Lambda (Xi dghter);#theta_{#Lambda};Number of #Lambda^0", 200, -10, 190);
 811      fListHistCascade->Add(fHistThetaLambdaXiMinus);
 812   }
 813
 814   if (!fHistThetaBachXiMinus) {
 815      fHistThetaBachXiMinus = new TH1F("fHistThetaBachXiMinus", "#theta of gen. Bach.;#theta_{Bach};Number of Bach.", 200, -10, 190);
 816      fListHistCascade->Add(fHistThetaBachXiMinus);
 817   }
 818
 819   if (!fHistThetaMesDghterXiMinus) {
 820      fHistThetaMesDghterXiMinus = new TH1F("fHistThetaMesDghterXiMinus", "#theta of gen. Meson #Lambda dghter;#theta_{MesDght};Number of Mes.", 200, -10, 190);
 821      fListHistCascade->Add(fHistThetaMesDghterXiMinus);
 822   }
 823
 824   if (!fHistThetaBarDghterXiMinus) {
 825      fHistThetaBarDghterXiMinus = new TH1F("fHistThetaBarDghterXiMinus", "#theta of gen. Baryon #Lambda dghter;#theta_{BarDght};Number of Bar.", 200, -10, 190);
 826      fListHistCascade->Add(fHistThetaBarDghterXiMinus);
 827   }
 828
 829                 // - Pt distribution (control plots)
 830
 831   if (!fHistPtBachXiMinus) {
 832      fHistPtBachXiMinus = new TH1F("fHistPtBachXiMinus", "p_{t} of gen. Bach.;pt_{Bach};Number of Bach.", 200, 0, 10);
 833      fListHistCascade->Add(fHistPtBachXiMinus);
 834   }
 835
 836   if (!fHistPtMesDghterXiMinus) {
 837      fHistPtMesDghterXiMinus = new TH1F("fHistPtMesDghterXiMinus", "p_{t} of gen. Meson #Lambda dghter;pt_{MesDght};Number of Mes.", 200, 0, 10);
 838      fListHistCascade->Add(fHistPtMesDghterXiMinus);
 839   }
 840
 841   if (!fHistPtBarDghterXiMinus) {
 842      fHistPtBarDghterXiMinus = new TH1F("fHistPtBarDghterXiMinus", "p_{t} of gen. Baryon #Lambda dghter;pt_{BarDght};Number of Bar.", 200, 0, 10);
 843      fListHistCascade->Add(fHistPtBarDghterXiMinus);
 844   }
 845
 846
 847
 848   //--------
 849   // II - Xi+
 850   // - Pseudo-Rapidity distribution
 851   if (!fHistEtaGenCascXiPlus) {
 852      fHistEtaGenCascXiPlus = new TH1F("fHistEtaGenCascXiPlus", "#eta of any gen. #Xi^{+};#eta;Number of Casc", 200, -10, 10);
 853      fListHistCascade->Add(fHistEtaGenCascXiPlus);
 854   }
 855   if (!f3dHistGenPtVsGenYvsCentXiPlusNat) {
 856      f3dHistGenPtVsGenYvsCentXiPlusNat = new TH3D("f3dHistGenPtVsGenYvsCentXiPlusNat", "MC P_{t} Vs MC Y of Gen #Xi^{+} ;Pt_{MC} (GeV/c); Y_{MC}", 100, ptBinLimits, 110, yBinLimits, 11, centBinLimits);
 857      fListHistCascade->Add(f3dHistGenPtVsGenYvsCentXiPlusNat);
 858   }
 859   if (!f3dHistGenPtVsGenYvsNtracksXiPlusNat) {
 860     f3dHistGenPtVsGenYvsNtracksXiPlusNat = new TH3D("f3dHistGenPtVsGenYvsNtracksXiPlusNat", "MC P_{t} Vs MC Y of Gen #Xi^{+} ;Pt_{MC} (GeV/c); Y_{MC}", 100, 0., 10., 110, -1.1, 1.1, 100, 0., 4000.);
 861     fListHistCascade->Add(f3dHistGenPtVsGenYvsNtracksXiPlusNat);
 862   }
 863   if (!f3dHistGenPtVsGenYvsCentXiPlusInj) {
 864      f3dHistGenPtVsGenYvsCentXiPlusInj = new TH3D("f3dHistGenPtVsGenYvsCentXiPlusInj", "MC P_{t} Vs MC Y of Gen #Xi^{+} ;Pt_{MC} (GeV/c); Y_{MC}", 100, ptBinLimits, 110, yBinLimits, 11, centBinLimits);
 865      fListHistCascade->Add(f3dHistGenPtVsGenYvsCentXiPlusInj);
 866   }
 867   if (!f3dHistGenPtVsGenYvsNtracksXiPlusInj) {
 868     f3dHistGenPtVsGenYvsNtracksXiPlusInj = new TH3D("f3dHistGenPtVsGenYvsNtracksXiPlusInj", "MC P_{t} Vs MC Y of Gen #Xi^{+} ;Pt_{MC} (GeV/c); Y_{MC}", 100, 0., 10., 110, -1.1, 1.1, 100, 0., 4000.);
 869     fListHistCascade->Add(f3dHistGenPtVsGenYvsNtracksXiPlusInj);
 870   }
 871   if (!f3dHistGenPtVsGenctauvsYXiPlusNat) {
 872      f3dHistGenPtVsGenctauvsYXiPlusNat = new TH3D("f3dHistGenPtVsGenctauvsYXiPlusNat", "MC P_{t} Vs MC ctau Vs Yof Gen #Xi^{+} ", 100, ptBinLimits, 111, ctauBinLimits, 110, yBinLimits);
 873      fListHistCascade->Add(f3dHistGenPtVsGenctauvsYXiPlusNat);
 874   }
 875   if (!f3dHistGenPtVsGenctauvsYXiPlusInj) {
 876      f3dHistGenPtVsGenctauvsYXiPlusInj = new TH3D("f3dHistGenPtVsGenctauvsYXiPlusInj", "MC P_{t} Vs MC ctau Vs Yof Gen #Xi^{+} ", 100, ptBinLimits, 111, ctauBinLimits, 110, yBinLimits);
 877      fListHistCascade->Add(f3dHistGenPtVsGenctauvsYXiPlusInj);
 878   }
 879
 880
 881
 882   // - Info at the generation level of multi-strange particle
 883
 884   if (!fHistThetaGenCascXiPlusNat) {
 885      fHistThetaGenCascXiPlusNat = new TH1F("fHistThetaGenCascXiPlusNat", "#theta of gen. #Xi^{+};#theta;Number of Casc.", 200, -10, 190);
 886      fListHistCascade->Add(fHistThetaGenCascXiPlusNat);
 887   }
 888   if (!fHistThetaGenCascXiPlusInj) {
 889      fHistThetaGenCascXiPlusInj = new TH1F("fHistThetaGenCascXiPlusInj", "#theta of inj. #Xi^{+};#theta;Number of Casc.", 200, -10, 190);
 890      fListHistCascade->Add(fHistThetaGenCascXiPlusInj);
 891   }
 892
 893
 894   if (!f2dHistGenPtVsGenYFdblXiPlus) {
 895      f2dHistGenPtVsGenYFdblXiPlus = new TH2D("f2dHistGenPtVsGenYFdblXiPlus", "MC P_{t} Vs MC Y of findable Gen #Xi^{+} ;Pt_{MC} (GeV/c); Y_{MC}", 100, 0., 10., 220, -1.1, 1.1);
 896      fListHistCascade->Add(f2dHistGenPtVsGenYFdblXiPlus);
 897   }
 898
 899                 // - Theta distribution the daughters (control plots)
 900
 901   if (!fHistThetaLambdaXiPlus) {
 902      fHistThetaLambdaXiPlus = new TH1F("fHistThetaLambdaXiPlus", "#theta of gen. #Lambda (Xi dghter);#theta_{#Lambda};Number of #Lambda", 200, -10, 190);
 903      fListHistCascade->Add(fHistThetaLambdaXiPlus);
 904   }
 905
 906   if (!fHistThetaBachXiPlus) {
 907      fHistThetaBachXiPlus = new TH1F("fHistThetaBachXiPlus", "#theta of gen. Bach.;#theta_{Bach};Number of Bach.", 200, -10, 190);
 908      fListHistCascade->Add(fHistThetaBachXiPlus);
 909   }
 910
 911   if (!fHistThetaMesDghterXiPlus) {
 912      fHistThetaMesDghterXiPlus = new TH1F("fHistThetaMesDghterXiPlus", "#theta of gen. Meson #Lambda dghter;#theta_{MesDght};Number of Mes.", 200, -10, 190);
 913      fListHistCascade->Add(fHistThetaMesDghterXiPlus);
 914   }
 915
 916   if (!fHistThetaBarDghterXiPlus) {
 917      fHistThetaBarDghterXiPlus = new TH1F("fHistThetaBarDghterXiPlus", "#theta of gen. Baryon #Lambda dghter;#theta_{BarDght};Number of Bar.", 200, -10, 190);
 918      fListHistCascade->Add(fHistThetaBarDghterXiPlus);
 919   }
 920
 921                 // - Pt distribution (control plots)
 922
 923   if (!fHistPtBachXiPlus) {
 924      fHistPtBachXiPlus = new TH1F("fHistPtBachXiPlus", "p_{t} of gen. Bach.;pt_{Bach};Number of Bach.", 200, 0, 10);
 925      fListHistCascade->Add(fHistPtBachXiPlus);
 926   }
 927
 928   if (!fHistPtMesDghterXiPlus) {
 929      fHistPtMesDghterXiPlus = new TH1F("fHistPtMesDghterXiPlus", "p_{t} of gen. Meson #Lambda dghter);pt_{MesDght};Number of Mes.", 200, 0, 10);
 930      fListHistCascade->Add(fHistPtMesDghterXiPlus);
 931   }
 932
 933   if (!fHistPtBarDghterXiPlus) {
 934      fHistPtBarDghterXiPlus = new TH1F("fHistPtBarDghterXiPlus", "p_{t} of gen. Baryon #Lambda dghter);pt_{BarDght};Number of Bar.", 200, 0, 10);
 935      fListHistCascade->Add(fHistPtBarDghterXiPlus);
 936   }
 937
 938
 939   //---------
 940   // III - Omega-
 941   // - Pseudo-Rapidity distribution
 942   if (!fHistEtaGenCascOmegaMinus) {
 943      fHistEtaGenCascOmegaMinus = new TH1F("fHistEtaGenCascOmegaMinus", "#eta of any gen. #Omega^{-};#eta;Number of Casc", 200, -10, 10);
 944      fListHistCascade->Add(fHistEtaGenCascOmegaMinus);
 945   }
 946  if (!f3dHistGenPtVsGenYvsCentOmegaMinusNat) {
 947      f3dHistGenPtVsGenYvsCentOmegaMinusNat = new TH3D("f3dHistGenPtVsGenYvsCentOmegaMinusNat", "MC P_{t} Vs MC Y of Gen #Omega^{-} ;Pt_{MC} (GeV/c); Y_{MC}", 100, ptBinLimits, 110, yBinLimits, 11, centBinLimits);
 948      fListHistCascade->Add(f3dHistGenPtVsGenYvsCentOmegaMinusNat);
 949   }
 950   if (!f3dHistGenPtVsGenYvsNtracksOmegaMinusNat) {
 951      f3dHistGenPtVsGenYvsNtracksOmegaMinusNat = new TH3D("f3dHistGenPtVsGenYvsNtracksOmegaMinusNat", "MC P_{t} Vs MC Y of Gen #Omega^{-} ;Pt_{MC} (GeV/c); Y_{MC}", 100, 0., 10., 110, -1.1, 1.1, 100, 0., 4000.);
 952      fListHistCascade->Add(f3dHistGenPtVsGenYvsNtracksOmegaMinusNat);
 953   }
 954   if (!f3dHistGenPtVsGenYvsCentOmegaMinusInj) {
 955      f3dHistGenPtVsGenYvsCentOmegaMinusInj = new TH3D("f3dHistGenPtVsGenYvsCentOmegaMinusInj", "MC P_{t} Vs MC Y of Gen #Omega^{-} ;Pt_{MC} (GeV/c); Y_{MC}", 100, ptBinLimits, 110, yBinLimits, 11, centBinLimits);
 956      fListHistCascade->Add(f3dHistGenPtVsGenYvsCentOmegaMinusInj);
 957   }
 958   if (!f3dHistGenPtVsGenYvsNtracksOmegaMinusInj) {
 959      f3dHistGenPtVsGenYvsNtracksOmegaMinusInj = new TH3D("f3dHistGenPtVsGenYvsNtracksOmegaMinusInj", "MC P_{t} Vs MC Y of Gen #Omega^{-} ;Pt_{MC} (GeV/c); Y_{MC}", 100, 0., 10., 110, -1.1, 1.1, 100, 0., 4000.);
 960      fListHistCascade->Add(f3dHistGenPtVsGenYvsNtracksOmegaMinusInj);
 961   }
 962   if (!f3dHistGenPtVsGenctauvsYOmegaMinusNat) {
 963      f3dHistGenPtVsGenctauvsYOmegaMinusNat = new TH3D("f3dHistGenPtVsGenctauvsYOmegaMinusNat", "MC P_{t} Vs MC ctau Vs Y of Gen #Omega^{-} ", 100, ptBinLimits, 111, ctauBinLimits, 110, yBinLimits);
 964      fListHistCascade->Add(f3dHistGenPtVsGenctauvsYOmegaMinusNat);
 965   }
 966   if (!f3dHistGenPtVsGenctauvsYOmegaMinusInj) {
 967      f3dHistGenPtVsGenctauvsYOmegaMinusInj = new TH3D("f3dHistGenPtVsGenctauvsYOmegaMinusInj", "MC P_{t} Vs MC ctau Vs Y of Gen #Omega^{-} ", 100, ptBinLimits, 111, ctauBinLimits, 110, yBinLimits);
 968      fListHistCascade->Add(f3dHistGenPtVsGenctauvsYOmegaMinusInj);
 969   }
 970
 971
 972
 973
 974
 975   // - Info at the generation level of multi-strange particle
 976
 977   if (!fHistThetaGenCascOmegaMinusNat) {
 978      fHistThetaGenCascOmegaMinusNat = new TH1F("fHistThetaGenCascOmegaMinusNat", "#theta of gen. #Omega^{-};#theta;Number of Casc.", 200, -10, 190);
 979      fListHistCascade->Add(fHistThetaGenCascOmegaMinusNat);
 980   }
 981   if (!fHistThetaGenCascOmegaMinusInj) {
 982      fHistThetaGenCascOmegaMinusInj = new TH1F("fHistThetaGenCascOmegaMinusInj", "#theta of inj. #Omega^{-};#theta;Number of Casc.", 200, -10, 190);
 983      fListHistCascade->Add(fHistThetaGenCascOmegaMinusInj);
 984   }
 985
 986
 987   if (!f2dHistGenPtVsGenYFdblOmegaMinus) {
 988      f2dHistGenPtVsGenYFdblOmegaMinus = new TH2D("f2dHistGenPtVsGenYFdblOmegaMinus", "MC P_{t} Vs MC Y of findable Gen #Omega^{-}; Pt_{MC} (GeV/c); Y_{MC}", 100, 0., 10., 220, -1.1, 1.1);
 989      fListHistCascade->Add(f2dHistGenPtVsGenYFdblOmegaMinus);
 990   }
 991
 992                 // - Theta distribution the daughters (control plots)
 993
 994   if (!fHistThetaLambdaOmegaMinus) {
 995      fHistThetaLambdaOmegaMinus = new TH1F("fHistThetaLambdaOmegaMinus", "#theta of gen. #Lambda (Omega dghter);#theta_{#Lambda};Number of #Lambda", 200, -10, 190);
 996      fListHistCascade->Add(fHistThetaLambdaOmegaMinus);
 997   }
 998
 999   if (!fHistThetaBachOmegaMinus) {
1000      fHistThetaBachOmegaMinus = new TH1F("fHistThetaBachOmegaMinus", "#theta of gen. Bach.;#theta_{Bach};Number of Bach.", 200, -10, 190);
1001      fListHistCascade->Add(fHistThetaBachOmegaMinus);
1002   }
1003
1004   if (!fHistThetaMesDghterOmegaMinus) {
1005      fHistThetaMesDghterOmegaMinus = new TH1F("fHistThetaMesDghterOmegaMinus", "#theta of gen. Meson #Lambda dghter;#theta_{MesDght};Number of Mes.", 200, -10, 190);
1006      fListHistCascade->Add(fHistThetaMesDghterOmegaMinus);
1007   }
1008
1009   if (!fHistThetaBarDghterOmegaMinus) {
1010      fHistThetaBarDghterOmegaMinus = new TH1F("fHistThetaBarDghterOmegaMinus", "#theta of gen. Baryon #Lambda dghter;#theta_{BarDght};Number of Bar.", 200, -10, 190);
1011      fListHistCascade->Add(fHistThetaBarDghterOmegaMinus);
1012   }
1013
1014                 // - Pt distribution (control plots)
1015
1016   if (!fHistPtBachOmegaMinus) {
1017      fHistPtBachOmegaMinus = new TH1F("fHistPtBachOmegaMinus", "p_{t} of gen. Bach.;pt_{Bach};Number of Bach.", 200, 0, 10);
1018      fListHistCascade->Add(fHistPtBachOmegaMinus);
1019   }
1020
1021   if (!fHistPtMesDghterOmegaMinus) {
1022      fHistPtMesDghterOmegaMinus = new TH1F("fHistPtMesDghterOmegaMinus", "p_{t} of gen. Meson #Lambda dghter);pt_{MesDght};Number of Mes.", 200, 0, 10);
1023      fListHistCascade->Add(fHistPtMesDghterOmegaMinus);
1024   }
1025
1026   if (!fHistPtBarDghterOmegaMinus) {
1027      fHistPtBarDghterOmegaMinus = new TH1F("fHistPtBarDghterOmegaMinus", "p_{t} of gen. Baryon #Lambda dghter);pt_{BarDght};Number of Bar.", 200, 0, 10);
1028      fListHistCascade->Add(fHistPtBarDghterOmegaMinus);
1029   }
1030
1031
1032   //---------
1033   // IV - Omega+
1034   // - Pseudo-Rapidity distribution
1035   if (!fHistEtaGenCascOmegaPlus) {
1036      fHistEtaGenCascOmegaPlus = new TH1F("fHistEtaGenCascOmegaPlus", "#eta of any gen. #Omega^{+};#eta;Number of Casc", 200, -10, 10);
1037      fListHistCascade->Add(fHistEtaGenCascOmegaPlus);
1038   }
1039   if (!f3dHistGenPtVsGenYvsCentOmegaPlusNat) {
1040      f3dHistGenPtVsGenYvsCentOmegaPlusNat = new TH3D("f3dHistGenPtVsGenYvsCentOmegaPlusNat", "MC P_{t} Vs MC Y of Gen #Omega^{+} ;Pt_{MC} (GeV/c); Y_{MC}", 100, ptBinLimits, 110, yBinLimits, 11, centBinLimits);
1041      fListHistCascade->Add(f3dHistGenPtVsGenYvsCentOmegaPlusNat);
1042   }
1043   if (!f3dHistGenPtVsGenYvsNtracksOmegaPlusNat) {
1044      f3dHistGenPtVsGenYvsNtracksOmegaPlusNat = new TH3D("f3dHistGenPtVsGenYvsNtracksOmegaPlusNat", "MC P_{t} Vs MC Y of Gen #Omega^{+} ;Pt_{MC} (GeV/c); Y_{MC}", 100, 0., 10., 110, -1.1, 1.1, 100, 0., 4000.);
1045      fListHistCascade->Add(f3dHistGenPtVsGenYvsNtracksOmegaPlusNat);
1046   }
1047   if (!f3dHistGenPtVsGenYvsCentOmegaPlusInj) {
1048      f3dHistGenPtVsGenYvsCentOmegaPlusInj = new TH3D("f3dHistGenPtVsGenYvsCentOmegaPlusInj", "MC P_{t} Vs MC Y of Gen #Omega^{+} ;Pt_{MC} (GeV/c); Y_{MC}", 100, ptBinLimits, 110, yBinLimits, 11, centBinLimits);
1049      fListHistCascade->Add(f3dHistGenPtVsGenYvsCentOmegaPlusInj);
1050   }
1051   if (!f3dHistGenPtVsGenYvsNtracksOmegaPlusInj) {
1052      f3dHistGenPtVsGenYvsNtracksOmegaPlusInj = new TH3D("f3dHistGenPtVsGenYvsNtracksOmegaPlusInj", "MC P_{t} Vs MC Y of Gen #Omega^{+} ;Pt_{MC} (GeV/c); Y_{MC}", 100, 0., 10., 110, -1.1, 1.1, 100, 0., 4000.);
1053      fListHistCascade->Add(f3dHistGenPtVsGenYvsNtracksOmegaPlusInj);
1054   }
1055   if (!f3dHistGenPtVsGenctauvsYOmegaPlusNat) {
1056      f3dHistGenPtVsGenctauvsYOmegaPlusNat = new TH3D("f3dHistGenPtVsGenctauvsYOmegaPlusNat", "MC P_{t} Vs MC ctau Vs Y of Gen #Omega^{+} ", 100, ptBinLimits, 111, ctauBinLimits, 110, yBinLimits);
1057      fListHistCascade->Add(f3dHistGenPtVsGenctauvsYOmegaPlusNat);
1058   }
1059   if (!f3dHistGenPtVsGenctauvsYOmegaPlusInj) {
1060      f3dHistGenPtVsGenctauvsYOmegaPlusInj = new TH3D("f3dHistGenPtVsGenctauvsYOmegaPlusInj", "MC P_{t} Vs MC ctau Vs Y of Gen #Omega^{+} ", 100, ptBinLimits, 111, ctauBinLimits, 110, yBinLimits);
1061      fListHistCascade->Add(f3dHistGenPtVsGenctauvsYOmegaPlusInj);
1062   }
1063
1064
1065   // - Info at the generation level of multi-strange particle
1066
1067   if (!fHistThetaGenCascOmegaPlusNat) {
1068      fHistThetaGenCascOmegaPlusNat = new TH1F("fHistThetaGenCascOmegaPlusNat", "#theta of gen. #Omega^{+};#theta;Number of Casc.", 200, -10, 190);
1069      fListHistCascade->Add(fHistThetaGenCascOmegaPlusNat);
1070   }
1071   if (!fHistThetaGenCascOmegaPlusInj) {
1072      fHistThetaGenCascOmegaPlusInj = new TH1F("fHistThetaGenCascOmegaPlusInj", "#theta of inj. #Omega^{+};#theta;Number of Casc.", 200, -10, 190);
1073      fListHistCascade->Add(fHistThetaGenCascOmegaPlusInj);
1074   }
1075
1076
1077   if (!f2dHistGenPtVsGenYFdblOmegaPlus) {
1078      f2dHistGenPtVsGenYFdblOmegaPlus = new TH2D("f2dHistGenPtVsGenYFdblOmegaPlus", "MC P_{t} Vs MC Y of findable Gen #Omega^{+}; Pt_{MC} (GeV/c); Y_{MC}", 100, 0., 10., 220, -1.1, 1.1);
1079      fListHistCascade->Add(f2dHistGenPtVsGenYFdblOmegaPlus);
1080   }
1081
1082
1083                 // - Theta distribution the daughters (control plots)
1084
1085   if (!fHistThetaLambdaOmegaPlus) {
1086      fHistThetaLambdaOmegaPlus = new TH1F("fHistThetaLambdaOmegaPlus", "#theta of gen. #Lambda (Omega dghter);#theta_{#Lambda};Number of #Lambda", 200, -10, 190);
1087      fListHistCascade->Add(fHistThetaLambdaOmegaPlus);
1088   }
1089
1090   if (!fHistThetaBachOmegaPlus) {
1091      fHistThetaBachOmegaPlus = new TH1F("fHistThetaBachOmegaPlus", "#theta of gen. Bach.;#theta_{Bach};Number of Bach.", 200, -10, 190);
1092      fListHistCascade->Add(fHistThetaBachOmegaPlus);
1093   }
1094
1095   if (!fHistThetaMesDghterOmegaPlus) {
1096      fHistThetaMesDghterOmegaPlus = new TH1F("fHistThetaMesDghterOmegaPlus", "#theta of gen. Meson #Lambda dghter;#theta_{MesDght};Number of Mes.", 200, -10, 190);
1097      fListHistCascade->Add(fHistThetaMesDghterOmegaPlus);
1098   }
1099
1100   if (!fHistThetaBarDghterOmegaPlus) {
1101      fHistThetaBarDghterOmegaPlus = new TH1F("fHistThetaBarDghterOmegaPlus", "#theta of gen. Baryon #Lambda dghter;#theta_{BarDght};Number of Bar.", 200, -10, 190);
1102      fListHistCascade->Add(fHistThetaBarDghterOmegaPlus);
1103   }
1104
1105                 // - Pt distribution (control plots)
1106
1107   if (!fHistPtBachOmegaPlus) {
1108      fHistPtBachOmegaPlus = new TH1F("fHistPtBachOmegaPlus", "p_{t} of gen. Bach.;pt_{Bach};Number of Bach.", 200, 0, 10);
1109      fListHistCascade->Add(fHistPtBachOmegaPlus);
1110   }
1111
1112   if (!fHistPtMesDghterOmegaPlus) {
1113      fHistPtMesDghterOmegaPlus = new TH1F("fHistPtMesDghterOmegaPlus", "p_{t} of gen. Meson #Lambda dghter);pt_{MesDght};Number of Mes.", 200, 0, 10);
1114      fListHistCascade->Add(fHistPtMesDghterOmegaPlus);
1115   }
1116
1117   if (!fHistPtBarDghterOmegaPlus) {
1118      fHistPtBarDghterOmegaPlus = new TH1F("fHistPtBarDghterOmegaPlus", "p_{t} of gen. Baryon #Lambda dghter);pt_{BarDght};Number of Bar.", 200, 0, 10);
1119      fListHistCascade->Add(fHistPtBarDghterOmegaPlus);
1120   }
1121
1122
1123 //--------------------------------------------------------------------------------
1124 // Part 2 - Any reconstructed cascades + reconstructed cascades associated with MC
1125
1126                 // - Effective mass histos for cascades candidates.
1127
1128   if (! fHistMassXiMinus) {
1129           fHistMassXiMinus = new TH1F("fHistMassXiMinus","#Xi^{-} candidates;M( #Lambda , #pi^{-} ) (GeV/c^{2});Counts", 400,1.2,2.0);
1130           fListHistCascade->Add(fHistMassXiMinus);
1131   }
1132
1133   if (! fHistMassXiPlus) {
1134           fHistMassXiPlus = new TH1F("fHistMassXiPlus","#Xi^{+} candidates;M( #bar{#Lambda}^{0} , #pi^{+} ) (GeV/c^{2});Counts",400,1.2,2.0);
1135           fListHistCascade->Add(fHistMassXiPlus);
1136   }
1137
1138   if (! fHistMassOmegaMinus) {
1139           fHistMassOmegaMinus = new TH1F("fHistMassOmegaMinus","#Omega^{-} candidates;M( #Lambda , K^{-} ) (GeV/c^{2});Counts", 500,1.5,2.5);
1140           fListHistCascade->Add(fHistMassOmegaMinus);
1141   }
1142
1143   if (! fHistMassOmegaPlus) {
1144           fHistMassOmegaPlus = new TH1F("fHistMassOmegaPlus","#Omega^{+} candidates;M( #bar{#Lambda}^{0} , K^{+} ) (GeV/c^{2});Counts", 500,1.5,2.5);
1145           fListHistCascade->Add(fHistMassOmegaPlus);
1146   }
1147
1148
1149
1150                 // - Effective mass histos with combined PID
1151
1152   if (! fHistMassWithCombPIDXiMinus) {
1153     fHistMassWithCombPIDXiMinus = new TH1F("fHistMassWithCombPIDXiMinus","#Xi^{-} candidates, with Bach. comb. PID;M( #Lambda , #pi^{-} ) (GeV/c^{2});Counts", 400,1.2,2.0);
1154     fListHistCascade->Add(fHistMassWithCombPIDXiMinus);
1155   }
1156
1157   if (! fHistMassWithCombPIDXiPlus) {
1158     fHistMassWithCombPIDXiPlus = new TH1F("fHistMassWithCombPIDXiPlus","#Xi^{+} candidates, with Bach. comb. PID;M( #bar{#Lambda}^{0} , #pi^{+} ) (GeV/c^{2});Counts",400,1.2,2.0);
1159     fListHistCascade->Add(fHistMassWithCombPIDXiPlus);
1160   }
1161
1162   if (! fHistMassWithCombPIDOmegaMinus) {
1163         fHistMassWithCombPIDOmegaMinus = new TH1F("fHistMassWithCombPIDOmegaMinus","#Omega^{-} candidates, with Bach. comb. PID;M( #Lambda , K^{-} ) (GeV/c^{2});Counts", 500,1.5,2.5);
1164     fListHistCascade->Add(fHistMassWithCombPIDOmegaMinus);
1165   }
1166
1167   if (! fHistMassWithCombPIDOmegaPlus) {
1168         fHistMassWithCombPIDOmegaPlus = new TH1F("fHistMassWithCombPIDOmegaPlus","#Omega^{+} candidates, with Bach. comb. PID;M( #bar{#Lambda}^{0} , K^{+} ) (GeV/c^{2});Counts", 500,1.5,2.5);
1169     fListHistCascade->Add(fHistMassWithCombPIDOmegaPlus);
1170   }
1171
1172                 // - PID Probability versus MC Pt(bachelor track)
1173   if(! f2dHistPIDprobaKaonVsMCPtBach ){
1174         f2dHistPIDprobaKaonVsMCPtBach  = new TH2F( "f2dHistPIDprobaKaonVsMCPtBach" , "Comb. PID proba to be K^{#pm} Vs MC Bach. Pt ; Pt_{MC}(Bach.) (GeV/c); Comb. PID Proba (Bach. = K^{#pm})", 100, 0.0, 5.0, 110, 0.0, 1.10 );
1175         fListHistCascade->Add(f2dHistPIDprobaKaonVsMCPtBach);
1176   }
1177
1178   if(! f2dHistPIDprobaPionVsMCPtBach ){
1179         f2dHistPIDprobaPionVsMCPtBach  = new TH2F( "f2dHistPIDprobaPionVsMCPtBach" , "Comb. PID proba to be #pi^{#pm} Vs MC Bach. Pt ; Pt_{MC}(Bach.) (GeV/c); Comb. PID Proba (Bach. = #pi^{#pm})", 100, 0.0, 5.0, 110, 0.0, 1.10 );
1180         fListHistCascade->Add(f2dHistPIDprobaPionVsMCPtBach);
1181   }
1182
1183
1184                 // - Effective mass histos with perfect MC PID on the bachelor
1185
1186   if (! fHistMassWithMcPIDXiMinus) {
1187     fHistMassWithMcPIDXiMinus = new TH1F("fHistMassWithMcPIDXiMinus","#Xi^{-} candidates, with Bach. MC PID;M( #Lambda , #pi^{-} ) (GeV/c^{2});Counts", 400,1.2,2.0);
1188     fListHistCascade->Add(fHistMassWithMcPIDXiMinus);
1189   }
1190
1191   if (! fHistMassWithMcPIDXiPlus) {
1192     fHistMassWithMcPIDXiPlus = new TH1F("fHistMassWithMcPIDXiPlus","#Xi^{+} candidates, with Bach. MC PID;M( #bar{#Lambda}^{0} , #pi^{+} ) (GeV/c^{2});Counts",400,1.2,2.0);
1193     fListHistCascade->Add(fHistMassWithMcPIDXiPlus);
1194   }
1195
1196   if (! fHistMassWithMcPIDOmegaMinus) {
1197         fHistMassWithMcPIDOmegaMinus = new TH1F("fHistMassWithMcPIDOmegaMinus","#Omega^{-} candidates, with Bach. MC PID;M( #Lambda , K^{-} ) (GeV/c^{2});Counts", 500,1.5,2.5);
1198     fListHistCascade->Add(fHistMassWithMcPIDOmegaMinus);
1199   }
1200
1201   if (! fHistMassWithMcPIDOmegaPlus) {
1202         fHistMassWithMcPIDOmegaPlus = new TH1F("fHistMassWithMcPIDOmegaPlus","#Omega^{+} candidates, with Bach. MC PID;M( #bar{#Lambda}^{0} , K^{+} ) (GeV/c^{2});Counts", 500,1.5,2.5);
1203     fListHistCascade->Add(fHistMassWithMcPIDOmegaPlus);
1204   }
1205
1206
1207                 // - Effective mass histos for cascades candidates ASSOCIATED with MC.
1208
1209   if (! fHistAsMCMassXiMinus) {
1210           fHistAsMCMassXiMinus = new TH1F("fHistAsMCMassXiMinus","#Xi^{-} candidates associated to MC;M( #Lambda , #pi^{-} ) (GeV/c^{2});Counts", 400,1.2,2.0);
1211           fListHistCascade->Add(fHistAsMCMassXiMinus);
1212   }
1213
1214   if (! fHistAsMCMassXiPlus) {
1215           fHistAsMCMassXiPlus = new TH1F("fHistAsMCMassXiPlus","#Xi^{+} candidates associated to MC;M( #bar{#Lambda}^{0} , #pi^{+} ) (GeV/c^{2});Counts",400,1.2,2.0);
1216           fListHistCascade->Add(fHistAsMCMassXiPlus);
1217   }
1218
1219   if (! fHistAsMCMassOmegaMinus) {
1220           fHistAsMCMassOmegaMinus = new TH1F("fHistAsMCMassOmegaMinus","#Omega^{-} candidates associated to MC;M( #Lambda , K^{-} ) (GeV/c^{2});Counts", 500,1.5,2.5);
1221           fListHistCascade->Add(fHistAsMCMassOmegaMinus);
1222   }
1223
1224   if (! fHistAsMCMassOmegaPlus) {
1225           fHistAsMCMassOmegaPlus = new TH1F("fHistAsMCMassOmegaPlus","#Omega^{+} candidates associated to MC;M( #bar{#Lambda}^{0} , K^{+} ) (GeV/c^{2});Counts", 500,1.5,2.5);
1226           fListHistCascade->Add(fHistAsMCMassOmegaPlus);
1227   }
1228
1229
1230                 // -  Generated Pt Vs generated Y of the cascade candidates associated with MC
1231                 //     + having the proper maximum proba of combined PID for the bachelor
1232
1233   if (!f2dHistAsMCandCombPIDGenPtVsGenYXiMinus) {
1234      f2dHistAsMCandCombPIDGenPtVsGenYXiMinus = new TH2F("f2dHistAsMCandCombPIDGenPtVsGenYXiMinus", "MC P_{t} Vs MC Y of #Xi^{-} (associated+Bach.PID); Pt_{MC} (GeV/c); Y_{MC}", 200, 0., 10., 220, -1.1, 1.1);
1235      fListHistCascade->Add(f2dHistAsMCandCombPIDGenPtVsGenYXiMinus);
1236   }
1237
1238   if (!f2dHistAsMCandCombPIDGenPtVsGenYXiPlus) {
1239      f2dHistAsMCandCombPIDGenPtVsGenYXiPlus = new TH2F("f2dHistAsMCandCombPIDGenPtVsGenYXiPlus", "MC P_{t} Vs MC Y of #Xi^{+} (associated+Bach.PID); Pt_{MC} (GeV/c); Y_{MC}", 100, 0., 10., 220, -1.1, 1.1);
1240      fListHistCascade->Add(f2dHistAsMCandCombPIDGenPtVsGenYXiPlus);
1241   }
1242
1243   if (!f2dHistAsMCandCombPIDGenPtVsGenYOmegaMinus) {
1244      f2dHistAsMCandCombPIDGenPtVsGenYOmegaMinus = new TH2F("f2dHistAsMCandCombPIDGenPtVsGenYOmegaMinus", "MC P_{t} Vs MC Y of #Omega^{-} (associated+Bach.PID); Pt_{MC} (GeV/c); Y_{MC}", 200, 0., 10., 220, -1.1, 1.1);
1245      fListHistCascade->Add(f2dHistAsMCandCombPIDGenPtVsGenYOmegaMinus);
1246   }
1247
1248   if (!f2dHistAsMCandCombPIDGenPtVsGenYOmegaPlus) {
1249      f2dHistAsMCandCombPIDGenPtVsGenYOmegaPlus = new TH2F("f2dHistAsMCandCombPIDGenPtVsGenYOmegaPlus", "MC P_{t} Vs MC Y of #Omega^{+} (associated+Bach.PID); Pt_{MC} (GeV/c); Y_{MC}", 200, 0., 10., 220, -1.1, 1.1);
1250      fListHistCascade->Add(f2dHistAsMCandCombPIDGenPtVsGenYOmegaPlus);
1251   }
1252
1253
1254                 // - Generated Pt Vs Generated Y, for the cascade candidates associated with MC
1255
1256   if (!f2dHistAsMCGenPtVsGenYXiMinus) {
1257         f2dHistAsMCGenPtVsGenYXiMinus = new TH2F("f2dHistAsMCGenPtVsGenYXiMinus", "MC P_{t} Vs MC Y of gen. #Xi^{-} (associated);Pt_{MC} (GeV/c); Rapidity, Y_{MC}",200, 0., 10., 220, -1.1, 1.1);
1258           fListHistCascade->Add(f2dHistAsMCGenPtVsGenYXiMinus );
1259   }
1260
1261   if (!f2dHistAsMCGenPtVsGenYXiPlus) {
1262           f2dHistAsMCGenPtVsGenYXiPlus = new TH2F("f2dHistAsMCGenPtVsGenYXiPlus", "MC P_{t} Vs MC Y of gen. #Xi^{+} (associated);Pt_{MC} (GeV/c); Rapidity, Y_{MC}",200, 0., 10., 220, -1.1, 1.1);
1263           fListHistCascade->Add(f2dHistAsMCGenPtVsGenYXiPlus );
1264   }
1265
1266   if (!f2dHistAsMCGenPtVsGenYOmegaMinus) {
1267           f2dHistAsMCGenPtVsGenYOmegaMinus = new TH2F("f2dHistAsMCGenPtVsGenYOmegaMinus", "MC P_{t} Vs MC Y of gen. #Omega^{-} (associated);Pt_{MC} (GeV/c); Rapidity, Y_{MC}",200, 0., 10., 220, -1.1, 1.1);
1268           fListHistCascade->Add(f2dHistAsMCGenPtVsGenYOmegaMinus );
1269   }
1270
1271   if (!f2dHistAsMCGenPtVsGenYOmegaPlus) {
1272           f2dHistAsMCGenPtVsGenYOmegaPlus = new TH2F("f2dHistAsMCGenPtVsGenYOmegaPlus", "MC P_{t} Vs MC Y of gen. #Omega^{+} (associated);Pt_{MC} (GeV/c); Rapidity, Y_{MC}",200, 0., 10., 220, -1.1, 1.1);
1273           fListHistCascade->Add(f2dHistAsMCGenPtVsGenYOmegaPlus );
1274   }
1275
1276
1277                 // - Generated Eta of the the cascade candidates associated with MC
1278   if (!fHistAsMCGenEtaXiMinus) {
1279           fHistAsMCGenEtaXiMinus = new TH1F("fHistAsMCGenEtaXiMinus", "#eta of gen. #Xi^{-} (associated);#eta;Number of Casc", 100, -5, 5);
1280           fListHistCascade->Add( fHistAsMCGenEtaXiMinus );
1281   }
1282
1283   if (!fHistAsMCGenEtaXiPlus) {
1284           fHistAsMCGenEtaXiPlus = new TH1F("fHistAsMCGenEtaXiPlus", "#eta of gen. #Xi^{+} (associated);#eta;Number of Casc", 100, -5, 5);
1285           fListHistCascade->Add( fHistAsMCGenEtaXiPlus );
1286   }
1287
1288   if (!fHistAsMCGenEtaOmegaMinus) {
1289           fHistAsMCGenEtaOmegaMinus = new TH1F("fHistAsMCGenEtaOmegaMinus", "#eta of gen. #Omega^{-} (associated);#eta;Number of Casc", 100, -5, 5);
1290           fListHistCascade->Add( fHistAsMCGenEtaOmegaMinus );
1291   }
1292
1293   if (!fHistAsMCGenEtaOmegaPlus) {
1294           fHistAsMCGenEtaOmegaPlus = new TH1F("fHistAsMCGenEtaOmegaPlus", "#eta of gen. #Omega^{+} (associated);#eta;Number of Casc", 100, -5, 5);
1295           fListHistCascade->Add( fHistAsMCGenEtaOmegaPlus );
1296   }
1297
1298
1299
1300                 // - Resolution in Pt as function of generated Pt
1301
1302   if(! f2dHistAsMCResPtXiMinus) {
1303           f2dHistAsMCResPtXiMinus = new TH2F( "f2dHistAsMCResPtXiMinus", "Resolution in Pt reconstruction for #Xi^{-}; Pt_{MC} (GeV/c); (Pt_{reco} - Pt_{MC}) / Pt_{MC}", 200, 0., 10., 200, -0.1, 0.1);
1304           fListHistCascade->Add(f2dHistAsMCResPtXiMinus);
1305   }
1306
1307   if(! f2dHistAsMCResPtXiPlus) {
1308           f2dHistAsMCResPtXiPlus = new TH2F( "f2dHistAsMCResPtXiPlus", "Resolution in Pt reconstruction for #Xi^{+}; Pt_{MC} (GeV/c); (Pt_{reco} - Pt_{MC}) / Pt_{MC}", 200, 0., 10., 200, -0.1, 0.1);
1309           fListHistCascade->Add(f2dHistAsMCResPtXiPlus);
1310   }
1311
1312   if(! f2dHistAsMCResPtOmegaMinus) {
1313           f2dHistAsMCResPtOmegaMinus = new TH2F( "f2dHistAsMCResPtOmegaMinus", "Resolution in Pt reconstruction for #Omega^{-}; Pt_{MC} (GeV/c); (Pt_{reco} - Pt_{MC}) / Pt_{MC}", 200, 0., 10., 200, -0.1, 0.1);
1314           fListHistCascade->Add(f2dHistAsMCResPtOmegaMinus);
1315   }
1316
1317   if(! f2dHistAsMCResPtOmegaPlus) {
1318           f2dHistAsMCResPtOmegaPlus = new TH2F( "f2dHistAsMCResPtOmegaPlus", "Resolution in Pt reconstruction for #Omega^{+}; Pt_{MC} (GeV/c); (Pt_{reco} - Pt_{MC}) / Pt_{MC}", 200, 0., 10., 200, -0.1, 0.1);
1319           fListHistCascade->Add(f2dHistAsMCResPtOmegaPlus);
1320   }
1321
1322                 // - Resolution in R(2D) as function of generated R
1323
1324   if(! f2dHistAsMCResRXiMinus) {
1325           f2dHistAsMCResRXiMinus = new TH2F( "f2dHistAsMCResRXiMinus", "Resolution in transv. position for #Xi^{-}; R_{MC} (cm); (R_{reco} - R_{MC}) / R_{MC}", 450, 0., 45.0, 240, -0.3, 0.3);
1326           fListHistCascade->Add(f2dHistAsMCResRXiMinus);
1327   }
1328
1329   if(! f2dHistAsMCResRXiPlus) {
1330           f2dHistAsMCResRXiPlus = new TH2F( "f2dHistAsMCResRXiPlus", "Resolution in transv. position for #Xi^{+}; R_{MC} (cm); (R_{reco} - R_{MC}) / R_{MC}", 450, 0., 45.0, 240, -0.3, 0.3);
1331           fListHistCascade->Add(f2dHistAsMCResRXiPlus);
1332   }
1333
1334   if(! f2dHistAsMCResROmegaMinus) {
1335           f2dHistAsMCResROmegaMinus = new TH2F( "f2dHistAsMCResROmegaMinus", "Resolution in transv. position for #Omega^{-}; R_{MC} (cm); (R_{reco} - R_{MC}) / R_{MC}", 450, 0., 45.0, 240, -0.3, 0.3);
1336           fListHistCascade->Add(f2dHistAsMCResROmegaMinus);
1337   }
1338
1339   if(! f2dHistAsMCResROmegaPlus) {
1340           f2dHistAsMCResROmegaPlus = new TH2F( "f2dHistAsMCResROmegaPlus", "Resolution in transv. position for #Omega^{+}; R_{MC} (cm); (R_{reco} - R_{MC}) / R_{MC}", 450, 0., 45.0, 240, -0.3, 0.3);
1341           fListHistCascade->Add(f2dHistAsMCResROmegaPlus);
1342   }
1343
1344                 // - Resolution in phi as function of generated Pt
1345
1346   if(! f2dHistAsMCResPhiXiMinus) {
1347           f2dHistAsMCResPhiXiMinus = new TH2F( "f2dHistAsMCResPhiXiMinus", "Resolution in #phi for #Xi^{-}; Pt_{MC} (GeV/c); #phi(MC) - #phi(reco)   (deg)", 200, 0., 10., 60, -30., 30.);
1348           fListHistCascade->Add(f2dHistAsMCResPhiXiMinus);
1349   }
1350
1351   if(! f2dHistAsMCResPhiXiPlus) {
1352           f2dHistAsMCResPhiXiPlus = new TH2F( "f2dHistAsMCResPhiXiPlus", "Resolution in #phi for #Xi^{+}; Pt_{MC} (GeV/c); #phi(MC) - #phi(reco)   (deg)", 200, 0., 10., 60, -30., 30.);
1353           fListHistCascade->Add(f2dHistAsMCResPhiXiPlus);
1354   }
1355
1356   if(! f2dHistAsMCResPhiOmegaMinus) {
1357           f2dHistAsMCResPhiOmegaMinus = new TH2F( "f2dHistAsMCResPhiOmegaMinus", "Resolution in #phi for #Omega^{-}; Pt_{MC} (GeV/c); #phi(MC) - #phi(reco)   (deg)", 200, 0., 10., 60, -30., 30.);
1358           fListHistCascade->Add(f2dHistAsMCResPhiOmegaMinus);
1359   }
1360
1361   if(! f2dHistAsMCResPhiOmegaPlus) {
1362           f2dHistAsMCResPhiOmegaPlus = new TH2F( "f2dHistAsMCResPhiOmegaPlus", "Resolution in #phi for #Omega^{+}; Pt_{MC} (GeV/c); #phi(MC) - #phi(reco)   (deg)", 200, 0., 10., 60, -30., 30.);
1363           fListHistCascade->Add(f2dHistAsMCResPhiOmegaPlus);
1364   }
1365
1366          //  - Correlation between proton (antiproton) daughter MC pt and Xi/Omega MC pt (to apply Geant/Fluka correction)
1367     if (!f2dHistAsMCptProtonMCptXiMinus) {
1368       f2dHistAsMCptProtonMCptXiMinus = new TH2F( "f2dHistAsMCptProtonMCptXiMinus", "Proton MC pt vs Xi- MC pt", 100, 0., 10., 100, 0., 10.);
1369       fListHistCascade->Add(f2dHistAsMCptProtonMCptXiMinus);
1370     }
1371     if (!f2dHistAsMCptAntiprotonMCptXiPlus) {
1372       f2dHistAsMCptAntiprotonMCptXiPlus = new TH2F( "f2dHistAsMCptAntiprotonMCptXiPlus", "Antiproton MC pt vs Xi+ MC pt", 100, 0., 10., 100, 0., 10.);
1373       fListHistCascade->Add(f2dHistAsMCptAntiprotonMCptXiPlus);
1374     }
1375     if (!f2dHistAsMCptProtonMCptOmegaMinus) {
1376       f2dHistAsMCptProtonMCptOmegaMinus = new TH2F( "f2dHistAsMCptProtonMCptOmegaMinus", "Proton MC pt vs Omega- MC pt", 100, 0., 10., 100, 0., 10.);
1377       fListHistCascade->Add(f2dHistAsMCptProtonMCptOmegaMinus);
1378     }
1379     if (!f2dHistAsMCptAntiprotonMCptOmegaPlus) {
1380       f2dHistAsMCptAntiprotonMCptOmegaPlus = new TH2F( "f2dHistAsMCptAntiprotonMCptOmegaPlus", "Antiproton MC pt vs Omega+ MC pt", 100, 0., 10., 100, 0., 10.);
1381       fListHistCascade->Add(f2dHistAsMCptAntiprotonMCptOmegaPlus);
1382     }
1383
1384     if (! fHistV0toXiCosineOfPointingAngle) {
1385         fHistV0toXiCosineOfPointingAngle = new TH1F("fHistV0toXiCosineOfPointingAngle", "Cos. of V0 Ptng Angl / Xi vtx ;Cos(V0 Point. Angl / Xi vtx); Counts", 200, 0.95, 1.0001);
1386         fListHistCascade->Add(fHistV0toXiCosineOfPointingAngle);
1387     }
1388
1389     if (! fHistV0CosineOfPointingAnglevsPtXi) {
1390         fHistV0CosineOfPointingAnglevsPtXi = new TH2F("fHistV0CosineOfPointingAnglevsPtXi", "Cos. of V0 Ptng Angl vs cascade Pt ;Cascade p_{T} (GeV/c); Cos(V0PA) wrt Xi vtx",100, 0., 10., 251, 0.95, 1.0002);
1391         fListHistCascade->Add(fHistV0CosineOfPointingAnglevsPtXi);
1392     }
1393     if (! fHistV0CosineOfPointingAnglevsPtOmega) {
1394         fHistV0CosineOfPointingAnglevsPtOmega = new TH2F("fHistV0CosineOfPointingAnglevsPtOmega", "Cos. of V0 Ptng Angl vs cascade Pt ;Cascade p_{T} (GeV/c); Cos(V0PA) wrt Omega vtx",100, 0., 10., 251, 0.95, 1.0002);
1395         fListHistCascade->Add(fHistV0CosineOfPointingAnglevsPtOmega);
1396     }
1397
1398     // - PID container
1399 if(! fCFContCascadePIDAsXiMinus)  {
1400   const Int_t  lNbSteps      =  7 ;
1401   const Int_t  lNbVariables  =  4 ;
1402
1403   //array for the number of bins in each dimension :
1404   Int_t lNbBinsPerVar[4] = {0};
1405   lNbBinsPerVar[0] = 100;
1406   lNbBinsPerVar[1] = 800;
1407   lNbBinsPerVar[2] = 22;
1408   lNbBinsPerVar[3] = 11;
1409
1410   fCFContCascadePIDAsXiMinus = new AliCFContainer("fCFContCascadePIDAsXiMinus","Pt_{cascade} Vs M_{#Xi^{-} candidates} Vs Y_{#Xi}", lNbSteps, lNbVariables, lNbBinsPerVar );
1411
1412   //setting the bin limits
1413   fCFContCascadePIDAsXiMinus->SetBinLimits(0,   0.0  ,  10.0 ); // Pt(Cascade)
1414   fCFContCascadePIDAsXiMinus->SetBinLimits(1,   1.2  ,   2.0 ); // Xi Effective mass
1415   fCFContCascadePIDAsXiMinus->SetBinLimits(2,  -1.1  ,   1.1 ); // Rapidity
1416   Double_t *lBinLim3  = new Double_t[ lNbBinsPerVar[3]+1 ];
1417   for(Int_t i=3; i< lNbBinsPerVar[3]+1;i++)   lBinLim3[i]  = (Double_t)(i-1)*10.;
1418   lBinLim3[0] = 0.0;
1419   lBinLim3[1] = 5.0;
1420   lBinLim3[2] = 10.0;
1421   fCFContCascadePIDAsXiMinus->SetBinLimits(3,  lBinLim3 );       // Centrality
1422
1423   // Setting the step title : one per PID case
1424   fCFContCascadePIDAsXiMinus->SetStepTitle(0, "No PID");
1425   fCFContCascadePIDAsXiMinus->SetStepTitle(1, "TPC PID / 4-#sigma cut on Bachelor track");
1426   fCFContCascadePIDAsXiMinus->SetStepTitle(2, "TPC PID / 4-#sigma cut on Bachelor+Baryon tracks");
1427   fCFContCascadePIDAsXiMinus->SetStepTitle(3, "TPC PID / 4-#sigma cut on Bachelor+Baryon+Meson tracks");
1428   fCFContCascadePIDAsXiMinus->SetStepTitle(4, "Comb. PID / Bachelor");
1429   fCFContCascadePIDAsXiMinus->SetStepTitle(5, "Comb. PID / Bachelor+Baryon");
1430   fCFContCascadePIDAsXiMinus->SetStepTitle(6, "Comb. PID / Bachelor+Baryon+Meson");
1431
1432   // Setting the variable title, per axis
1433   fCFContCascadePIDAsXiMinus->SetVarTitle(0, "Pt_{cascade} (GeV/c)");
1434   fCFContCascadePIDAsXiMinus->SetVarTitle(1, "M( #Lambda , #pi^{-} ) (GeV/c^{2})");
1435   fCFContCascadePIDAsXiMinus->SetVarTitle(2, "Y_{#Xi}");
1436   fCFContCascadePIDAsXiMinus->SetVarTitle(3, "Centrality");
1437
1438   fListHistCascade->Add(fCFContCascadePIDAsXiMinus);
1439
1440 }
1441
1442 if(! fCFContCascadePIDAsXiPlus)  {
1443   const Int_t  lNbSteps      =  7 ;
1444   const Int_t  lNbVariables  =  4 ;
1445
1446   //array for the number of bins in each dimension :
1447   Int_t lNbBinsPerVar[4] = {0};
1448   lNbBinsPerVar[0] = 100;
1449   lNbBinsPerVar[1] = 800;
1450   lNbBinsPerVar[2] = 22;
1451   lNbBinsPerVar[3] = 11;
1452
1453   fCFContCascadePIDAsXiPlus = new AliCFContainer("fCFContCascadePIDAsXiPlus","Pt_{cascade} Vs M_{#Xi^{+} candidates} Vs Y_{#Xi}", lNbSteps, lNbVariables, lNbBinsPerVar );
1454
1455
1456   //setting the bin limits (valid  for v4-18-10-AN)
1457   fCFContCascadePIDAsXiPlus->SetBinLimits(0,   0.0  ,  10.0 );  // Pt(Cascade)
1458   fCFContCascadePIDAsXiPlus->SetBinLimits(1,   1.2  ,   2.0 );  // Xi Effective mass
1459   fCFContCascadePIDAsXiPlus->SetBinLimits(2,  -1.1  ,   1.1 );  // Rapidity
1460   Double_t *lBinLim3  = new Double_t[ lNbBinsPerVar[3]+1 ];
1461   for(Int_t i=3; i< lNbBinsPerVar[3]+1;i++)   lBinLim3[i]  = (Double_t)(i-1)*10.;
1462   lBinLim3[0] = 0.0;
1463   lBinLim3[1] = 5.0;
1464   lBinLim3[2] = 10.0;
1465   fCFContCascadePIDAsXiPlus->SetBinLimits(3,lBinLim3);     // Centrality
1466
1467   // Setting the step title : one per PID case
1468   fCFContCascadePIDAsXiPlus->SetStepTitle(0, "No PID");
1469   fCFContCascadePIDAsXiPlus->SetStepTitle(1, "TPC PID / 4-#sigma cut on Bachelor track");
1470   fCFContCascadePIDAsXiPlus->SetStepTitle(2, "TPC PID / 4-#sigma cut on Bachelor+Baryon tracks");
1471   fCFContCascadePIDAsXiPlus->SetStepTitle(3, "TPC PID / 4-#sigma cut on Bachelor+Baryon+Meson tracks");
1472   fCFContCascadePIDAsXiPlus->SetStepTitle(4, "Comb. PID / Bachelor");
1473   fCFContCascadePIDAsXiPlus->SetStepTitle(5, "Comb. PID / Bachelor+Baryon");
1474   fCFContCascadePIDAsXiPlus->SetStepTitle(6, "Comb. PID / Bachelor+Baryon+Meson");
1475
1476   // Setting the variable title, per axis
1477   fCFContCascadePIDAsXiPlus->SetVarTitle(0, "Pt_{cascade} (GeV/c)");
1478   fCFContCascadePIDAsXiPlus->SetVarTitle(1, "M( #Lambda , #pi^{+} ) (GeV/c^{2})");
1479   fCFContCascadePIDAsXiPlus->SetVarTitle(2, "Y_{#Xi}");
1480   fCFContCascadePIDAsXiPlus->SetVarTitle(3, "Centrality");
1481
1482   fListHistCascade->Add(fCFContCascadePIDAsXiPlus);
1483
1484 }
1485
1486
1487 if(! fCFContCascadePIDAsOmegaMinus)  {
1488   const Int_t  lNbSteps      =  7 ;
1489   const Int_t  lNbVariables  =  4 ;
1490
1491   //array for the number of bins in each dimension :
1492   Int_t lNbBinsPerVar[4] = {0};
1493   lNbBinsPerVar[0] = 100;
1494   lNbBinsPerVar[1] = 1000;
1495   lNbBinsPerVar[2] = 22;
1496   lNbBinsPerVar[3] = 11;
1497
1498   fCFContCascadePIDAsOmegaMinus = new AliCFContainer("fCFContCascadePIDAsOmegaMinus","Pt_{cascade} Vs M_{#Omega^{-} candidates} Vs Y_{#Omega}", lNbSteps, lNbVariables, lNbBinsPerVar );
1499
1500
1501   //setting the bin limits
1502   fCFContCascadePIDAsOmegaMinus->SetBinLimits(0,   0.0  ,  10.0 );      // Pt(Cascade)
1503   fCFContCascadePIDAsOmegaMinus->SetBinLimits(1,   1.5  ,   2.5 );      // Omega Effective mass
1504   fCFContCascadePIDAsOmegaMinus->SetBinLimits(2,  -1.1  ,   1.1 );      // Rapidity
1505   Double_t *lBinLim3  = new Double_t[ lNbBinsPerVar[3]+1 ];
1506   for(Int_t i=3; i< lNbBinsPerVar[3]+1;i++)   lBinLim3[i]  = (Double_t)(i-1)*10.;
1507   lBinLim3[0] = 0.0;
1508   lBinLim3[1] = 5.0;
1509   lBinLim3[2] = 10.0;
1510   fCFContCascadePIDAsOmegaMinus->SetBinLimits(3,lBinLim3);     // Centrality
1511
1512   // Setting the step title : one per PID case
1513   fCFContCascadePIDAsOmegaMinus->SetStepTitle(0, "No PID");
1514   fCFContCascadePIDAsOmegaMinus->SetStepTitle(1, "TPC PID / 4-#sigma cut on Bachelor track");
1515   fCFContCascadePIDAsOmegaMinus->SetStepTitle(2, "TPC PID / 4-#sigma cut on Bachelor+Baryon tracks");
1516   fCFContCascadePIDAsOmegaMinus->SetStepTitle(3, "TPC PID / 4-#sigma cut on Bachelor+Baryon+Meson tracks");
1517   fCFContCascadePIDAsOmegaMinus->SetStepTitle(4, "Comb. PID / Bachelor");
1518   fCFContCascadePIDAsOmegaMinus->SetStepTitle(5, "Comb. PID / Bachelor+Baryon");
1519   fCFContCascadePIDAsOmegaMinus->SetStepTitle(6, "Comb. PID / Bachelor+Baryon+Meson");
1520
1521   // Setting the variable title, per axis
1522   fCFContCascadePIDAsOmegaMinus->SetVarTitle(0, "Pt_{cascade} (GeV/c)");
1523   fCFContCascadePIDAsOmegaMinus->SetVarTitle(1, "M( #Lambda , K^{-} ) (GeV/c^{2})");
1524   fCFContCascadePIDAsOmegaMinus->SetVarTitle(2, "Y_{#Omega}");
1525   fCFContCascadePIDAsOmegaMinus->SetVarTitle(3, "Centrality");
1526
1527   fListHistCascade->Add(fCFContCascadePIDAsOmegaMinus);
1528
1529 }
1530
1531 if(! fCFContCascadePIDAsOmegaPlus)  {
1532   const Int_t  lNbSteps      =  7 ;
1533   const Int_t  lNbVariables  =  4 ;
1534
1535   //array for the number of bins in each dimension :
1536   Int_t lNbBinsPerVar[4]= {0};
1537   lNbBinsPerVar[0] = 100;
1538   lNbBinsPerVar[1] = 1000;
1539   lNbBinsPerVar[2] = 22;
1540   lNbBinsPerVar[3] = 11;
1541
1542   fCFContCascadePIDAsOmegaPlus = new AliCFContainer("fCFContCascadePIDAsOmegaPlus","Pt_{cascade} Vs M_{#Omega^{+} candidates} Vs Y_{#Omega}", lNbSteps, lNbVariables, lNbBinsPerVar );
1543
1544
1545   //setting the bin limits
1546   fCFContCascadePIDAsOmegaPlus->SetBinLimits(0,   0.0  ,  10.0 );       // Pt(Cascade)
1547   fCFContCascadePIDAsOmegaPlus->SetBinLimits(1,   1.5  ,   2.5 );       // Omega Effective mass
1548   fCFContCascadePIDAsOmegaPlus->SetBinLimits(2,  -1.1  ,   1.1 );       // Rapidity
1549   Double_t *lBinLim3  = new Double_t[ lNbBinsPerVar[3]+1 ];
1550   for(Int_t i=3; i< lNbBinsPerVar[3]+1;i++)   lBinLim3[i]  = (Double_t)(i-1)*10.;
1551   lBinLim3[0] = 0.0;
1552   lBinLim3[1] = 5.0;
1553   lBinLim3[2] = 10.0;
1554   fCFContCascadePIDAsOmegaPlus->SetBinLimits(3,lBinLim3);     // Centrality
1555
1556   // Setting the step title : one per PID case
1557   fCFContCascadePIDAsOmegaPlus->SetStepTitle(0, "No PID");
1558   fCFContCascadePIDAsOmegaPlus->SetStepTitle(1, "TPC PID / 4-#sigma cut on Bachelor track");
1559   fCFContCascadePIDAsOmegaPlus->SetStepTitle(2, "TPC PID / 4-#sigma cut on Bachelor+Baryon tracks");
1560   fCFContCascadePIDAsOmegaPlus->SetStepTitle(3, "TPC PID / 4-#sigma cut on Bachelor+Baryon+Meson tracks");
1561   fCFContCascadePIDAsOmegaPlus->SetStepTitle(4, "Comb. PID / Bachelor");
1562   fCFContCascadePIDAsOmegaPlus->SetStepTitle(5, "Comb. PID / Bachelor+Baryon");
1563   fCFContCascadePIDAsOmegaPlus->SetStepTitle(6, "Comb. PID / Bachelor+Baryon+Meson");
1564
1565   // Setting the variable title, per axis
1566   fCFContCascadePIDAsOmegaPlus->SetVarTitle(0, "Pt_{cascade} (GeV/c)");
1567   fCFContCascadePIDAsOmegaPlus->SetVarTitle(1, "M( #Lambda , K^{+} ) (GeV/c^{2})");
1568   fCFContCascadePIDAsOmegaPlus->SetVarTitle(2, "Y_{#Omega}");
1569   fCFContCascadePIDAsOmegaPlus->SetVarTitle(3, "Centrality");
1570
1571   fListHistCascade->Add(fCFContCascadePIDAsOmegaPlus);
1572
1573 }
1574
1575 // Part 3 : Towards the optimisation of topological selections -------
1576 if(! fCFContAsCascadeCuts){
1577
1578         // Container meant to store all the relevant distributions corresponding to the cut variables.
1579         //          - NB overflow/underflow of variables on which we want to cut later should be 0!!!
1580
1581   const Int_t  lNbSteps      =  4 ;
1582   const Int_t  lNbVariables  =  22 ;
1583
1584   //array for the number of bins in each dimension :
1585   Int_t lNbBinsPerVar[lNbVariables] = {0};
1586   lNbBinsPerVar[0]  = 100;
1587   lNbBinsPerVar[1]  = 126;
1588   lNbBinsPerVar[2]  = 24;
1589   lNbBinsPerVar[3]  = 220;
1590   lNbBinsPerVar[4]  = 30;
1591   lNbBinsPerVar[5]  = 50;
1592
1593   lNbBinsPerVar[6]  = 101;
1594   lNbBinsPerVar[7]  = 102;
1595   lNbBinsPerVar[8]  = 101;
1596   lNbBinsPerVar[9]  = 26;
1597   lNbBinsPerVar[10] = 26;
1598
1599   lNbBinsPerVar[11] = 150; // 1-MeV/c2 bins
1600   lNbBinsPerVar[12] = 120;
1601
1602   lNbBinsPerVar[13] = 100;
1603
1604   lNbBinsPerVar[14] = 110; // 0.02 unit of y per bin
1605   lNbBinsPerVar[15] = 110;
1606
1607   lNbBinsPerVar[16] = 2;
1608
1609   lNbBinsPerVar[17] = 11;
1610   lNbBinsPerVar[18] = 100;
1611
1612   lNbBinsPerVar[19] = 112; // Proper time of cascade
1613   lNbBinsPerVar[20] = 112; // Proper time of V0
1614   lNbBinsPerVar[21] = 112; // Distance V0-Xi in the transverse plane
1615
1616   fCFContAsCascadeCuts = new AliCFContainer("fCFContAsCascadeCuts","Cut Container for Asso. Cascades", lNbSteps, lNbVariables, lNbBinsPerVar );
1617
1618   //0
1619   fCFContAsCascadeCuts->SetBinLimits(0, 0., 2.);                 // DcaXiDaughters : 0.0 to 2.0
1620   //1
1621    Double_t *lBinLim1  = new Double_t[ lNbBinsPerVar[1]+1 ];
1622    lBinLim1[0] = 0.0;
1623    lBinLim1[1] = 0.03;
1624    for(Int_t i=2; i< lNbBinsPerVar[1];i++)   lBinLim1[i]  = (Double_t)0.03   + (5.  - 0.03 )/(lNbBinsPerVar[1]-2)  * (Double_t)(i-1) ;
1625    lBinLim1[ lNbBinsPerVar[1]  ] = 100.0;
1626    fCFContAsCascadeCuts -> SetBinLimits(1,  lBinLim1 );
1627   delete [] lBinLim1;                                            // DcaBachToPrimVertexXi : 0.0 to 0.5
1628   //2
1629   fCFContAsCascadeCuts->SetBinLimits(2, .9988, 1.);              // XiCosineOfPointingAngle : 0.99 to 1.0
1630   //3
1631   fCFContAsCascadeCuts -> SetBinLimits(3, 0., 110. );            // XiRadius : 0.0 to 110.0
1632
1633   //4
1634   fCFContAsCascadeCuts->SetBinLimits(4, 1.1, 1.13);              // InvMassLambdaAsCascDghter
1635   //5
1636   fCFContAsCascadeCuts->SetBinLimits(5, 0., 2.);                 // DcaV0DaughtersXi : 0.0 to 2.0
1637   //6
1638   fCFContAsCascadeCuts->SetBinLimits(6, 0.98, 1.0002);           // V0CosineOfPointingAngleXi : 0.98 to 1.0
1639   //7
1640   Double_t *lBinLim7  = new Double_t[ lNbBinsPerVar[7]+1 ];
1641   for(Int_t i=0; i< lNbBinsPerVar[7]-1;i++)   lBinLim7[i]  = (Double_t)0.0   + (100.  - 0.0 )/(lNbBinsPerVar[7]-2)  * (Double_t)i ;
1642   lBinLim7[ lNbBinsPerVar[7]-1] = 200.0;
1643   lBinLim7[ lNbBinsPerVar[7]] = 1000.0;
1644   fCFContAsCascadeCuts -> SetBinLimits(7,  lBinLim7 );
1645   delete [] lBinLim7;                                             // V0RadiusXi : 0.0 to 1000.0
1646   //8
1647   Double_t *lBinLim8  = new Double_t[ lNbBinsPerVar[8]+1 ];
1648   for(Int_t i=0; i< lNbBinsPerVar[8];i++)   lBinLim8[i]  = (Double_t)0.0   + (0.4  - 0.0 )/(lNbBinsPerVar[8]-1)  * (Double_t)i ;
1649   lBinLim8[ lNbBinsPerVar[8]  ] = 100.0;
1650   fCFContAsCascadeCuts -> SetBinLimits(8,  lBinLim8 );
1651   delete [] lBinLim8;                                             // DcaV0ToPrimVertexXi : 0. to 0.4
1652   //9
1653   Double_t *lBinLim9  = new Double_t[ lNbBinsPerVar[9]+1 ];
1654   for(Int_t i=0; i< lNbBinsPerVar[9];i++)   lBinLim9[i]  = (Double_t)0.0   + (0.25  - 0.0 )/(lNbBinsPerVar[9]-1)  * (Double_t)i ;
1655   lBinLim9[ lNbBinsPerVar[9]  ] = 100.0;
1656   fCFContAsCascadeCuts -> SetBinLimits(9,  lBinLim9 );
1657   delete [] lBinLim9;                                             // DcaPosToPrimVertexXi : 0.0 to 0.25
1658   //10
1659   Double_t *lBinLim10  = new Double_t[ lNbBinsPerVar[10]+1 ];
1660   for(Int_t i=0; i< lNbBinsPerVar[10];i++)   lBinLim10[i]  = (Double_t)0.0   + (0.25  - 0.0 )/(lNbBinsPerVar[10]-1)  * (Double_t)i ;
1661   lBinLim10[ lNbBinsPerVar[10]  ] = 100.0;
1662   fCFContAsCascadeCuts -> SetBinLimits(10,  lBinLim10 );
1663   delete [] lBinLim10;                                            // DcaPosToPrimVertexXi : 0.0 to 0.25
1664   //11
1665   fCFContAsCascadeCuts->SetBinLimits(11, 1.25, 1.40);             // InvMassXi
1666   fCFContAsCascadeCuts->SetBinLimits(12, 1.62, 1.74);             // InvMassOmega
1667   fCFContAsCascadeCuts->SetBinLimits(13, 0.0, 10.0);              // XiTransvMom
1668   fCFContAsCascadeCuts->SetBinLimits(14, -1.1, 1.1);              // Y(Xi)
1669   fCFContAsCascadeCuts->SetBinLimits(15, -1.1, 1.1);              // Y(Omega)
1670   fCFContAsCascadeCuts->SetBinLimits(16, 0.0, 2.0);               // Natural or injected
1671   Double_t *lBinLim17  = new Double_t[ lNbBinsPerVar[17]+1 ];
1672   for(Int_t i=3; i< lNbBinsPerVar[17]+1;i++)   lBinLim17[i]  = (Double_t)(i-1)*10.;
1673   lBinLim17[0] = 0.0;
1674   lBinLim17[1] = 5.0;
1675   lBinLim17[2] = 10.0;
1676   fCFContAsCascadeCuts -> SetBinLimits(17,  lBinLim17 );       // Centrality
1677   delete [] lBinLim17;
1678   fCFContAsCascadeCuts->SetBinLimits(18, 0.0, 6000.0);         // ESD track multiplicity
1679   Double_t *lBinLim19  = new Double_t[ lNbBinsPerVar[19]+1 ];
1680   for(Int_t i=0; i< lNbBinsPerVar[19];i++)   lBinLim19[i]  = (Double_t)-1.   + (110.  + 1.0 )/(lNbBinsPerVar[19]-1)  * (Double_t)i ;
1681   lBinLim19[ lNbBinsPerVar[19]  ] = 2000.0;
1682   fCFContAsCascadeCuts->SetBinLimits(19, lBinLim19);           // Proper time cascade
1683   fCFContAsCascadeCuts->SetBinLimits(20, lBinLim19);           // Proper time V0
1684   fCFContAsCascadeCuts->SetBinLimits(21, lBinLim19);           // Distance V0-Xi in the transverse plane
1685   delete [] lBinLim19;
1686   // Setting the number of steps : one for each cascade species (Xi-, Xi+ and Omega-, Omega+)
1687   fCFContAsCascadeCuts->SetStepTitle(0, "#Xi^{-} candidates associated to MC");
1688   fCFContAsCascadeCuts->SetStepTitle(1, "#bar{#Xi}^{+} candidates associated to MC");
1689   fCFContAsCascadeCuts->SetStepTitle(2, "#Omega^{-} candidates associated to MC");
1690   fCFContAsCascadeCuts->SetStepTitle(3, "#bar{#Omega}^{+} candidates associated to MC");
1691
1692   // Setting the variable title, per axis
1693   fCFContAsCascadeCuts->SetVarTitle(0,  "DCA(XiDaughters) (cm)");
1694   fCFContAsCascadeCuts->SetVarTitle(1,  "DCA(Bach/PrimVertex) (cm)");
1695   fCFContAsCascadeCuts->SetVarTitle(2,  "cos(Xi pointing angle)");
1696   fCFContAsCascadeCuts->SetVarTitle(3,  "R_{2d}(Xi decay) (cm)");
1697   fCFContAsCascadeCuts->SetVarTitle(4,  "M_{#Lambda}(As Casc Dghter) (GeV/c^{2})");
1698   fCFContAsCascadeCuts->SetVarTitle(5,  "DCA(V0 Daughters Xi) (cm)");
1699
1700   fCFContAsCascadeCuts->SetVarTitle(6,  "cos(V0 pointing Angle) in Casc");
1701   fCFContAsCascadeCuts->SetVarTitle(7,  "R_{2d}(V0 decay) (cm)");
1702   fCFContAsCascadeCuts->SetVarTitle(8,  "DCA(V0/PrimVertex) (cm)");
1703   fCFContAsCascadeCuts->SetVarTitle(9,  "DCA(Pos/PrimVertex) (cm)");
1704   fCFContAsCascadeCuts->SetVarTitle(10, "DCA(Neg/PrimVertex) (cm)");
1705
1706   fCFContAsCascadeCuts->SetVarTitle(11, "Inv. Mass(Xi) (GeV/c^{2})");
1707   fCFContAsCascadeCuts->SetVarTitle(12, "Inv. Mass(Omega) (GeV/c^{2})");
1708
1709   fCFContAsCascadeCuts->SetVarTitle(13, "Pt_{MC}(Casc.) (GeV/c)");
1710
1711   fCFContAsCascadeCuts->SetVarTitle(14, "Y_{MC}(Xi)");
1712   fCFContAsCascadeCuts->SetVarTitle(15, "Y_{MC}(Omega)");
1713
1714   fCFContAsCascadeCuts->SetVarTitle(16, "Injected or natural particles");
1715
1716   fCFContAsCascadeCuts->SetVarTitle(17, "Centrality");
1717   fCFContAsCascadeCuts->SetVarTitle(18, "ESD track multiplicity");
1718
1719   fCFContAsCascadeCuts->SetVarTitle(19, "Proper time cascade");
1720   fCFContAsCascadeCuts->SetVarTitle(20, "Proper time V0");
1721   fCFContAsCascadeCuts->SetVarTitle(21, "Distance V0-Xi in the transverse plane");
1722
1723   fListHistCascade->Add(fCFContAsCascadeCuts);
1724 }
1725
1726   fV0Ampl = new TH1F("fV0Ampl","",500,0.,30000);
1727   fListHistCascade->Add(fV0Ampl);
1728
1729
1730 PostData(1, fListHistCascade);
1731 }// end CreateOutputObjects
1732
1733
1734 //________________________________________________________________________
1735 void AliAnalysisTaskCheckPerformanceCascadePbPb::UserExec(Option_t *) {
1736
1737   // Main loop
1738   // Called for each event
1739
1740         AliESDEvent *lESDevent = 0x0;
1741         AliAODEvent *lAODevent = 0x0;
1742         AliMCEvent  *lMCevent  = 0x0;
1743         AliStack    *lMCstack  = 0x0;
1744
1745         TClonesArray *arrayMC = 0;
1746
1747         if (!fPIDResponse) {
1748           AliError("Cannot get pid response");
1749           return;
1750         }
1751
1752
1753   // Connect to the InputEvent
1754   // After these lines, we should have an ESD/AOD event
1755
1756         if (fAnalysisType == "ESD") {
1757
1758           lESDevent = dynamic_cast<AliESDEvent*>( InputEvent() );
1759           if (!lESDevent) {
1760             Printf("ERROR: lESDevent not available \n");
1761             cout << "Name of the file with pb :" <<  CurrentFileName() << endl;
1762             return;
1763           }
1764           lMCevent = MCEvent();
1765           if (!lMCevent) {
1766             Printf("ERROR: Could not retrieve MC event \n");
1767             cout << "Name of the file with pb :" <<  CurrentFileName() << endl;
1768             return;
1769           }
1770
1771           lMCstack = lMCevent->Stack();
1772           if (!lMCstack) {
1773             Printf("ERROR: Could not retrieve MC stack \n");
1774             cout << "Name of the file with pb :" <<  CurrentFileName() << endl;
1775             return;
1776           }
1777
1778           if (fkRerunV0CascVertexers) { // relaunch V0 and Cascade vertexer
1779            /* lESDevent->ResetCascades();
1780             lESDevent->ResetV0s();
1781
1782             AliV0vertexer lV0vtxer;
1783             AliCascadeVertexer lCascVtxer;
1784
1785             lV0vtxer.SetCuts(fV0Sels);
1786             lCascVtxer.SetCuts(fCascSels);
1787
1788             lV0vtxer.Tracks2V0vertices(lESDevent);
1789             lCascVtxer.V0sTracks2CascadeVertices(lESDevent);*/
1790           }
1791
1792         } else if (fAnalysisType == "AOD") {
1793
1794           lAODevent = dynamic_cast<AliAODEvent*>( InputEvent() );
1795           if (!lAODevent) {
1796             Printf("ERROR: lAODevent not available \n");
1797             cout << "Name of the file with pb :" <<  CurrentFileName() << endl;
1798             return;
1799           }
1800
1801           arrayMC = (TClonesArray*) lAODevent->GetList()->FindObject(AliAODMCParticle::StdBranchName());
1802           if (!arrayMC) AliFatal("Error: MC particles branch not found!\n");
1803
1804         } else {
1805
1806           Printf("Analysis type (ESD or AOD) not specified \n");
1807           return;
1808
1809         }
1810
1811
1812   //------------------------------------------------
1813   // 1 - Preparing the general info about of the event = prim. Vtx + magnetic field (ESD)
1814
1815         Double_t lBestPrimaryVtxPos[3]   = {-100.0, -100.0, -100.0};
1816         Double_t lMagneticField    = -10.;
1817         Int_t ncascades = -1;
1818         // Primary tracks from ESD/AOD
1819         Float_t lPrimaryTrackMultiplicity= -1.;
1820
1821 //        Int_t    lSPDTrackletsMultiplicity = -1;
1822
1823         AliCentrality* centrality = 0x0;
1824         AliESDVZERO* esdV0 = 0x0;
1825         AliAODVZERO* aodV0 = 0x0;
1826         Float_t multV0A = 0.;
1827         Float_t multV0C = 0.;
1828
1829         Int_t   nNumberOfMCPrimaries     = -1;
1830
1831         Int_t nTrackMultiplicity = (InputEvent())->GetNumberOfTracks();
1832
1833         if (fAnalysisType == "ESD" ) {
1834
1835         // Prim vertex
1836           const AliESDVertex *lPrimaryTrackingVtx = lESDevent->GetPrimaryVertexTracks();  // get the vtx stored in ESD found with tracks
1837           const AliESDVertex *lPrimarySPDVtx      = lESDevent->GetPrimaryVertexSPD();     // get the vtx stored in ESD found with SPD tracklets
1838
1839           const AliESDVertex *lPrimaryBestVtx     = lESDevent->GetPrimaryVertex();
1840           // get the best primary vertex available for the event
1841           // As done in AliCascadeVertexer, we keep the one which is the best one available.
1842           // between : Tracking vertex > SPD vertex > TPC vertex > default SPD vertex
1843           lPrimaryBestVtx->GetXYZ( lBestPrimaryVtxPos );
1844
1845           // FIXME : quality selection regarding pile-up rejection, not in the Cascade task
1846           if (fkRejectEventPileUp) {
1847             if (lESDevent->IsPileupFromSPD() ) {// minContributors=3, minZdist=0.8, nSigmaZdist=3., nSigmaDiamXY=2., nSigmaDiamZ=5.  -> see http://alisoft.cern.ch/viewvc/trunk/STEER/AliESDEvent.h?root=AliRoot&r1=41914&r2=42199&pathrev=42199
1848               AliWarning("Pb / Event tagged as pile-up by SPD... return !");
1849               PostData(1, fListHistCascade);
1850               return;
1851             }
1852           }
1853           // FIXME : remove TPC-only primary vertex : retain only events with tracking + SPD vertex
1854           if (fkQualityCutNoTPConlyPrimVtx) {
1855             if (!lPrimarySPDVtx->GetStatus() && !lPrimaryTrackingVtx->GetStatus() ){
1856               AliWarning("Pb / No SPD prim. vertex nor prim. Tracking vertex ... return !");
1857               PostData(1, fListHistCascade);
1858               return;
1859             }
1860           }
1861
1862           lMagneticField = lESDevent->GetMagneticField();
1863
1864           ncascades = lESDevent->GetNumberOfCascades();
1865
1866           lPrimaryTrackMultiplicity = fESDtrackCuts->CountAcceptedTracks(lESDevent);
1867
1868 //          const AliMultiplicity *lAliMult = lESDevent->GetMultiplicity();
1869 //          lSPDTrackletsMultiplicity       = lAliMult->GetNumberOfTracklets();
1870
1871           //        nNumberOfMCPrimaries       = lMCstack->GetNprimary();
1872           nNumberOfMCPrimaries = lMCstack->GetNtrack(); // MN: this stack->GetNtrack(); has to be used to check primaries because in HIJING decay products of D and B mesons are also primaries and produced in HIJING during transport, however this is not the number of primaries!
1873           // Centrality determination
1874           esdV0 = lESDevent->GetVZEROData();
1875           multV0A = esdV0->GetMTotV0A();
1876           multV0C = esdV0->GetMTotV0C();
1877
1878           centrality = lESDevent->GetCentrality();
1879
1880         } else if (fAnalysisType == "AOD") {
1881
1882           const AliAODVertex *lPrimaryBestAODVtx = lAODevent->GetPrimaryVertex();
1883           if (!lPrimaryBestAODVtx) {
1884             AliWarning("No prim. vertex in AOD... return!");
1885             PostData(1, fListHistCascade);
1886             return;
1887           }
1888
1889           // get the best primary vertex available for the event GetVertex(0)
1890           // This one will be used for next calculations (DCA essentially)
1891           lPrimaryBestAODVtx->GetXYZ( lBestPrimaryVtxPos );
1892
1893           lMagneticField = lAODevent->GetMagneticField();
1894
1895           ncascades = lAODevent->GetNumberOfCascades();
1896
1897           lPrimaryTrackMultiplicity = 0.;
1898           for (Int_t itrack = 0; itrack<nTrackMultiplicity; itrack++) {
1899             AliAODTrack* track = lAODevent->GetTrack(itrack);
1900             if (track->TestFilterBit(AliAODTrack::kTrkGlobalNoDCA)) lPrimaryTrackMultiplicity++;
1901           }
1902
1903 //          lSPDTrackletsMultiplicity = lAODevent->GetTracklets()->GetNumberOfTracklets();
1904
1905           nNumberOfMCPrimaries = arrayMC->GetEntries();
1906
1907
1908           centrality = lAODevent->GetCentrality();
1909           aodV0 = lAODevent->GetVZEROData();
1910           multV0A=aodV0->GetMTotV0A();
1911           multV0C=aodV0->GetMTotV0C();
1912
1913         }
1914
1915         if (fkQualityCutZprimVtxPos) { // Should be already in the centrality selection
1916           if (TMath::Abs(lBestPrimaryVtxPos[2]) > fVtxRange ) {
1917             AliWarning("Pb / | Z position of Best Prim Vtx | > 10.0 cm ... return !");
1918             PostData(1, fListHistCascade);
1919             return;
1920           }
1921         }
1922
1923
1924       //  Printf("Centrality percentile V0M for this event %f)\n",  centrality->GetCentralityPercentile("V0M"));
1925
1926
1927         Float_t lcentrality = centrality->GetCentralityPercentile(fCentrEstimator.Data());
1928         if (lcentrality<fCentrLowLim||lcentrality>=fCentrUpLim) {
1929           PostData(1, fListHistCascade);
1930           return;
1931         }
1932
1933         if (nNumberOfMCPrimaries < 1) {
1934           PostData(1, fListHistCascade);
1935           return; // should be useless because we require vertex and centrality selection
1936         }
1937
1938         fV0Ampl->Fill(multV0A+multV0C);
1939
1940         // Best vertex distribution
1941         fHistBestVtxX->Fill(lBestPrimaryVtxPos[0]);
1942         fHistBestVtxY->Fill(lBestPrimaryVtxPos[1]);
1943         fHistBestVtxZ->Fill(lBestPrimaryVtxPos[2]);
1944
1945         fHistEvtsInCentralityBinsvsNtracks->Fill(lcentrality,lPrimaryTrackMultiplicity);
1946
1947         fHistMCTrackMultiplicity->Fill( nNumberOfMCPrimaries );  // MN: neutral particles included and also not physical ones
1948
1949         //cout << "Name of the accessed file :" <<  CurrentFileName() << endl;
1950
1951         Int_t   nMCPrimariesInAcceptance   =  0;
1952 //_____________________________________________________________________________
1953 // Part 1 - Loop over the MC primaries
1954
1955     for (Int_t iCurrentLabelStack = 0; iCurrentLabelStack < nNumberOfMCPrimaries; iCurrentLabelStack++) {
1956
1957       if (fAnalysisType == "ESD") {
1958
1959         TParticle* lCurrentParticle = 0x0;
1960                    lCurrentParticle = lMCstack->Particle( iCurrentLabelStack );
1961         if (!lCurrentParticle) {
1962           Printf("MC Primary loop %d - MC TParticle pointer to current stack particle = 0x0 ! Skip ...\n", iCurrentLabelStack );
1963           continue;
1964         }
1965
1966         if (!lMCstack->IsPhysicalPrimary(iCurrentLabelStack)) continue;
1967
1968
1969         if ( TMath::Abs( lCurrentParticle->Eta() ) > 0.8 ) continue;
1970
1971       } else if (fAnalysisType == "AOD") {
1972
1973         AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iCurrentLabelStack);
1974         if(TMath::Abs(partMC->Eta())>0.8) continue;
1975
1976       }
1977
1978       nMCPrimariesInAcceptance++;
1979
1980    }
1981
1982    f2dHistRecoMultVsMCMult->Fill( lPrimaryTrackMultiplicity, nMCPrimariesInAcceptance );  // MN: neutral are included
1983
1984
1985    // For proton
1986
1987 /*   for (Int_t iCurrentLabelStack = 0; iCurrentLabelStack < nNumberOfMCPrimaries; iCurrentLabelStack++) {
1988
1989         TParticle* lCurrentParticle = lMCstack->Particle( iCurrentLabelStack );
1990         if(!lCurrentParticle){
1991                 Printf("Proton loop %d - MC TParticle pointer to current stack particle = 0x0 ! Skip ...\n", iCurrentLabelStack );
1992                 continue;
1993
1994         }
1995
1996         if( lCurrentParticle->GetPdgCode() == 2212 )
1997                 fHistEtaGenProton->Fill( lCurrentParticle->Eta() );
1998
1999         if( lCurrentParticle->GetPdgCode() == -2212 )
2000                 fHistEtaGenAntiProton->Fill( lCurrentParticle->Eta() );
2001       }// end loop over primary proton
2002 */
2003
2004
2005
2006 //_____________________________________________________________________________
2007 // Part 2 - Loop over the different types of GENERATED cascades (Xi-+, Omega-+)
2008
2009         // - Initialisation of useful local variables
2010
2011         Int_t lPdgCodeCasc            = 0;
2012         Int_t lPdgCodeBach            = 0;
2013         Int_t lPdgCodeLambda          = 0;
2014         Int_t lPdgCodeDghtMesV0       = 0;
2015         Int_t lPdgCodeDghtBarV0       = 0;
2016
2017
2018         TH1F *lHistEtaGenCasc         = 0;
2019         TH3D *l3dHistGenPtVsGenYvsCentNat = 0;
2020         TH3D *l3dHistGenPtVsGenYvsNtracksNat = 0;
2021         TH3D *l3dHistGenPtVsGenYvsCentInj = 0;
2022         TH3D *l3dHistGenPtVsGenYvsNtracksInj = 0;
2023         TH3D *l3dHistGenPtVsGenctauvsYNat = 0;
2024         TH3D *l3dHistGenPtVsGenctauvsYInj = 0;
2025
2026         TH1F *lHistThetaGenCascNat    = 0;
2027         TH1F *lHistThetaGenCascInj    = 0;
2028         TH2D *l2dHistGenPtVsGenYFdbl  = 0;
2029         TH1F *lHistThetaLambda        = 0;
2030         TH1F *lHistThetaBach          = 0;
2031         TH1F *lHistThetaBarDghter     = 0;
2032         TH1F *lHistThetaMesDghter     = 0;
2033         TH1F *lHistPtBach             = 0;
2034         TH1F *lHistPtBarDghter        = 0;
2035         TH1F *lHistPtMesDghter        = 0;
2036         Int_t ncascperev = 0;
2037         Int_t ncascperevtot = 0;
2038
2039         Int_t endOfHijingEvent = -1;
2040
2041 for (Int_t iCascType = 1; iCascType < 5; iCascType++) {
2042   ncascperev = 0;
2043   ncascperevtot = 0;
2044
2045   endOfHijingEvent = -1;
2046
2047   switch (iCascType) {
2048     case 1: // Xi-
2049          lPdgCodeCasc       =   3312;  //Xi-
2050          lPdgCodeBach       =   -211;  //Pi-
2051          lPdgCodeLambda     =   3122;  //Lambda0
2052          lPdgCodeDghtMesV0  =   -211;  //Pi-
2053          lPdgCodeDghtBarV0  =   2212;  //Proton
2054
2055                 // any Xi-
2056          lHistEtaGenCasc        = fHistEtaGenCascXiMinus;
2057          l3dHistGenPtVsGenYvsCentNat = f3dHistGenPtVsGenYvsCentXiMinusNat;
2058          l3dHistGenPtVsGenYvsNtracksNat = f3dHistGenPtVsGenYvsNtracksXiMinusNat;
2059          l3dHistGenPtVsGenYvsCentInj = f3dHistGenPtVsGenYvsCentXiMinusInj;
2060          l3dHistGenPtVsGenYvsNtracksInj = f3dHistGenPtVsGenYvsNtracksXiMinusInj;
2061          l3dHistGenPtVsGenctauvsYNat = f3dHistGenPtVsGenctauvsYXiMinusNat;
2062          l3dHistGenPtVsGenctauvsYInj = f3dHistGenPtVsGenctauvsYXiMinusInj;
2063
2064                 // cascades generated within acceptance (cut in pt + theta)
2065          lHistThetaGenCascNat      = fHistThetaGenCascXiMinusNat;
2066          lHistThetaGenCascInj      = fHistThetaGenCascXiMinusInj;
2067          l2dHistGenPtVsGenYFdbl = f2dHistGenPtVsGenYFdblXiMinus;
2068          lHistThetaLambda       = fHistThetaLambdaXiMinus;
2069          lHistThetaBach         = fHistThetaBachXiMinus;
2070          lHistThetaBarDghter    = fHistThetaBarDghterXiMinus;
2071          lHistThetaMesDghter    = fHistThetaMesDghterXiMinus;
2072          lHistPtBach            = fHistPtBachXiMinus;
2073          lHistPtBarDghter       = fHistPtBarDghterXiMinus;
2074          lHistPtMesDghter       = fHistPtMesDghterXiMinus;
2075         break;
2076
2077     case 2: // Xi+
2078          lPdgCodeCasc        =  -3312;  //Xi+
2079          lPdgCodeBach        =    211;  //Pi+
2080          lPdgCodeLambda      =  -3122;  //AntiLambda0
2081          lPdgCodeDghtMesV0   =    211;  //Pi+
2082          lPdgCodeDghtBarV0   =  -2212;  //AntiProton
2083
2084                 // any Xi+
2085          lHistEtaGenCasc        = fHistEtaGenCascXiPlus;
2086          l3dHistGenPtVsGenYvsCentNat = f3dHistGenPtVsGenYvsCentXiPlusNat;
2087          l3dHistGenPtVsGenYvsNtracksNat = f3dHistGenPtVsGenYvsNtracksXiPlusNat;
2088          l3dHistGenPtVsGenYvsCentInj = f3dHistGenPtVsGenYvsCentXiPlusInj;
2089          l3dHistGenPtVsGenYvsNtracksInj = f3dHistGenPtVsGenYvsNtracksXiPlusInj;
2090          l3dHistGenPtVsGenctauvsYNat = f3dHistGenPtVsGenctauvsYXiPlusNat;
2091          l3dHistGenPtVsGenctauvsYInj = f3dHistGenPtVsGenctauvsYXiPlusInj;
2092
2093                 // cascades generated within acceptance (cut in pt + theta)
2094          lHistThetaGenCascNat      = fHistThetaGenCascXiPlusNat;
2095          lHistThetaGenCascInj      = fHistThetaGenCascXiPlusInj;
2096          l2dHistGenPtVsGenYFdbl = f2dHistGenPtVsGenYFdblXiPlus;
2097          lHistThetaLambda       = fHistThetaLambdaXiPlus;
2098          lHistThetaBach         = fHistThetaBachXiPlus;
2099          lHistThetaBarDghter    = fHistThetaBarDghterXiPlus;
2100          lHistThetaMesDghter    = fHistThetaMesDghterXiPlus;
2101          lHistPtBach            = fHistPtBachXiPlus;
2102          lHistPtBarDghter       = fHistPtBarDghterXiPlus;
2103          lHistPtMesDghter       = fHistPtMesDghterXiPlus;
2104         break;
2105
2106     case 3: // Omega-
2107          lPdgCodeCasc       =   3334;  //Omega-
2108          lPdgCodeBach       =   -321;  //K-
2109          lPdgCodeLambda     =   3122;  //Lambda0
2110          lPdgCodeDghtMesV0  =   -211;  //Pi-
2111          lPdgCodeDghtBarV0  =   2212;  //Proton
2112
2113                 // any Omega-
2114          lHistEtaGenCasc        = fHistEtaGenCascOmegaMinus;
2115          l3dHistGenPtVsGenYvsCentNat = f3dHistGenPtVsGenYvsCentOmegaMinusNat;
2116          l3dHistGenPtVsGenYvsNtracksNat = f3dHistGenPtVsGenYvsNtracksOmegaMinusNat;
2117          l3dHistGenPtVsGenYvsCentInj = f3dHistGenPtVsGenYvsCentOmegaMinusInj;
2118          l3dHistGenPtVsGenYvsNtracksInj = f3dHistGenPtVsGenYvsNtracksOmegaMinusInj;
2119          l3dHistGenPtVsGenctauvsYNat = f3dHistGenPtVsGenctauvsYOmegaMinusNat;
2120          l3dHistGenPtVsGenctauvsYInj = f3dHistGenPtVsGenctauvsYOmegaMinusInj;
2121
2122                 // cascades generated within acceptance (cut in pt + theta)
2123          lHistThetaGenCascNat      = fHistThetaGenCascOmegaMinusNat;
2124          lHistThetaGenCascInj      = fHistThetaGenCascOmegaMinusInj;
2125          l2dHistGenPtVsGenYFdbl = f2dHistGenPtVsGenYFdblOmegaMinus;
2126          lHistThetaLambda       = fHistThetaLambdaOmegaMinus;
2127          lHistThetaBach         = fHistThetaBachOmegaMinus;
2128          lHistThetaBarDghter    = fHistThetaBarDghterOmegaMinus;
2129          lHistThetaMesDghter    = fHistThetaMesDghterOmegaMinus;
2130          lHistPtBach            = fHistPtBachOmegaMinus;
2131          lHistPtBarDghter       = fHistPtBarDghterOmegaMinus;
2132          lHistPtMesDghter       = fHistPtMesDghterOmegaMinus;
2133         break;
2134
2135     case 4:  // Omega+
2136          lPdgCodeCasc       =  -3334;  //Omega+
2137          lPdgCodeBach       =    321;  //K+
2138          lPdgCodeLambda     =  -3122;  //AntiLambda0
2139          lPdgCodeDghtMesV0  =    211;  //Pi+
2140          lPdgCodeDghtBarV0  =  -2212;  //AntiProton
2141
2142                 // any Omega+
2143          lHistEtaGenCasc        = fHistEtaGenCascOmegaPlus;
2144          l3dHistGenPtVsGenYvsCentNat = f3dHistGenPtVsGenYvsCentOmegaPlusNat;
2145          l3dHistGenPtVsGenYvsNtracksNat = f3dHistGenPtVsGenYvsNtracksOmegaPlusNat;
2146          l3dHistGenPtVsGenYvsCentInj = f3dHistGenPtVsGenYvsCentOmegaPlusInj;
2147          l3dHistGenPtVsGenYvsNtracksInj = f3dHistGenPtVsGenYvsNtracksOmegaPlusInj;
2148          l3dHistGenPtVsGenctauvsYNat = f3dHistGenPtVsGenctauvsYOmegaPlusNat;
2149          l3dHistGenPtVsGenctauvsYInj = f3dHistGenPtVsGenctauvsYOmegaPlusInj;
2150
2151
2152                 // cascades generated within acceptance (cut in pt + theta)
2153          lHistThetaGenCascNat      = fHistThetaGenCascOmegaPlusNat;
2154          lHistThetaGenCascInj      = fHistThetaGenCascOmegaPlusInj;
2155          l2dHistGenPtVsGenYFdbl = f2dHistGenPtVsGenYFdblOmegaPlus;
2156          lHistThetaLambda       = fHistThetaLambdaOmegaPlus;
2157          lHistThetaBach         = fHistThetaBachOmegaPlus;
2158          lHistThetaBarDghter    = fHistThetaBarDghterOmegaPlus;
2159          lHistThetaMesDghter    = fHistThetaMesDghterOmegaPlus;
2160          lHistPtBach            = fHistPtBachOmegaPlus;
2161          lHistPtBarDghter       = fHistPtBarDghterOmegaPlus;
2162          lHistPtMesDghter       = fHistPtMesDghterOmegaPlus;
2163         break;
2164
2165   }// end switch cascade
2166
2167
2168   for (Int_t iCurrentLabelStack = 0; iCurrentLabelStack < nNumberOfMCPrimaries; iCurrentLabelStack++) {
2169
2170     Bool_t kIsNaturalPart = kFALSE;
2171     Double_t partEnergy = 0.;
2172     Double_t partPz     = 0.;
2173     Double_t partEta    = 0.;
2174     Double_t partTheta  = 0.;
2175     Double_t partP      = 0.;
2176     Double_t partPt     = 0.;
2177     Double_t partVx     = 0.;
2178     Double_t partVy     = 0.;
2179     Double_t partVz     = 0.;
2180     Double_t bacVx      = 0.;
2181     Double_t bacVy      = 0.;
2182     Double_t bacVz      = 0.;
2183     Double_t partMass   = 0.;
2184
2185     if ( fAnalysisType == "ESD" ) {
2186       TParticle* lCurrentParticle = 0x0;
2187                  lCurrentParticle = lMCstack->Particle( iCurrentLabelStack );
2188       if (!lCurrentParticle) {
2189         Printf("Cascade loop %d - MC TParticle pointer to current stack particle = 0x0 ! Skip ...\n", iCurrentLabelStack );
2190         continue;
2191       }
2192       if (!lMCstack->IsPhysicalPrimary(iCurrentLabelStack)) continue;
2193       if (lCurrentParticle->GetPdgCode() == lPdgCodeCasc) {  // Here !
2194         //cout << "Xi- within loop " << iCurrentLabelStack << "/ " << nNumberOfMCPrimaries << endl;
2195
2196         // -  Xi level ... _____________________________________________________________
2197         TParticle* xiMC = 0x0;
2198                    xiMC = lCurrentParticle;
2199         if (!xiMC) {
2200           Printf("MC TParticle pointer to Cascade = 0x0 ! Skip ...");
2201           continue;
2202         }
2203
2204         kIsNaturalPart = lMCevent->IsFromBGEvent(iCurrentLabelStack);
2205         if (!kIsNaturalPart) {
2206           if (!(lCurrentParticle->GetFirstMother()<0)) kIsNaturalPart = kTRUE; // because there are primaries (ALICE definition) not produced in the collision
2207         }
2208         partEnergy = xiMC->Energy();
2209         partPz = xiMC->Pz();
2210         partEta = xiMC->Eta();
2211         partPt = xiMC->Pt();
2212         partP = xiMC->P();
2213         partTheta = xiMC->Theta();
2214         partMass = xiMC->GetMass();
2215         partVx = xiMC->Vx();
2216         partVy = xiMC->Vy();
2217         partVz = xiMC->Vz();
2218
2219         if (xiMC->GetDaughter(0)>=0) {
2220           TParticle *mcBach = lMCstack->Particle(xiMC->GetDaughter(0));
2221           if (mcBach) {
2222             bacVx  = mcBach->Vx();
2223             bacVy  = mcBach->Vy();
2224             bacVz  = mcBach->Vz();
2225           }
2226         }
2227       } else continue;
2228     } else if ( fAnalysisType == "AOD" ) {
2229
2230       AliAODMCParticle *lCurrentParticleaod = (AliAODMCParticle*) arrayMC->At(iCurrentLabelStack);
2231       if (!lCurrentParticleaod) {
2232         Printf("Cascade loop %d - MC TParticle pointer to current stack particle = 0x0 ! Skip ...\n", iCurrentLabelStack );
2233         continue;
2234       }
2235
2236       if (endOfHijingEvent==-1) {  // For injected MC: determine where HIJING event ends
2237         if ((lCurrentParticleaod->GetStatus()==21)||
2238             ( (lCurrentParticleaod->GetPdgCode() == 443) &&
2239               (lCurrentParticleaod->GetMother() == -1)   &&
2240               (lCurrentParticleaod->GetDaughter(0) ==  (iCurrentLabelStack+1)))) {
2241
2242
2243           endOfHijingEvent = iCurrentLabelStack;
2244
2245         }
2246       }
2247
2248       if (!lCurrentParticleaod->IsPhysicalPrimary()) continue;
2249       if (!(lCurrentParticleaod->PdgCode() == lPdgCodeCasc)) continue;
2250
2251
2252       kIsNaturalPart = kTRUE;
2253       if ((iCurrentLabelStack>=endOfHijingEvent)&&(endOfHijingEvent!=-1)&&(lCurrentParticleaod->GetMother()<0)) kIsNaturalPart = kFALSE;
2254
2255       partEnergy = lCurrentParticleaod->E();
2256       partPz = lCurrentParticleaod->Pz();
2257       partEta = lCurrentParticleaod->Eta();
2258       partP = lCurrentParticleaod->P();
2259       partPt = lCurrentParticleaod->Pt();
2260       partTheta = lCurrentParticleaod->Theta();
2261       partMass = lCurrentParticleaod->M(); // FIXME not sure this works, seems not implemented
2262       partVx = lCurrentParticleaod->Xv();
2263       partVy = lCurrentParticleaod->Yv();
2264       partVz = lCurrentParticleaod->Zv();
2265
2266       if (lCurrentParticleaod->GetDaughter(0)>=0) {
2267         AliAODMCParticle *mcBach = (AliAODMCParticle*) arrayMC->At(lCurrentParticleaod->GetDaughter(0));
2268         if (mcBach) {
2269           bacVx  = mcBach->Xv();
2270           bacVy  = mcBach->Yv();
2271           bacVz  = mcBach->Zv();
2272         }
2273       }
2274
2275     }
2276
2277     ncascperevtot++;
2278         // Fill the first histos : = any generated Xi, not necessarily within the acceptance
2279     Double_t lRapXiMC = 0.5*TMath::Log((partEnergy + partPz) / (partEnergy - partPz +1.e-13));
2280         // Calculate proper time
2281     Double_t lctau = TMath::Sqrt((partVx-bacVx)*(partVx-bacVx)+(partVy-bacVy)*(partVy-bacVy)+(partVz-bacVz)*(partVz-bacVz));
2282     if (partP!=0.)    lctau = lctau*partMass/partP;
2283     else lctau = -1.;
2284
2285     Double_t lRadToDeg = 180.0/TMath::Pi();
2286
2287     // Fill the first histos : = any generated Xi, not necessarily within the acceptance
2288     lHistEtaGenCasc     ->Fill( partEta );
2289
2290     if (kIsNaturalPart) {
2291       l3dHistGenPtVsGenYvsCentNat   ->Fill( partPt, lRapXiMC, lcentrality );
2292       l3dHistGenPtVsGenYvsNtracksNat->Fill( partPt, lRapXiMC, lPrimaryTrackMultiplicity );
2293       l3dHistGenPtVsGenctauvsYNat->Fill( partPt, lctau, lRapXiMC );
2294       lHistThetaGenCascNat       ->Fill( lRadToDeg * partTheta );
2295     } else {
2296       l3dHistGenPtVsGenYvsCentInj   ->Fill( partPt, lRapXiMC, lcentrality );
2297       l3dHistGenPtVsGenYvsNtracksInj->Fill( partPt, lRapXiMC, lPrimaryTrackMultiplicity );
2298       l3dHistGenPtVsGenctauvsYInj->Fill( partPt, lctau, lRapXiMC );
2299       lHistThetaGenCascInj       ->Fill( lRadToDeg * partTheta );
2300     }
2301
2302     // Check the emission of particle stays within the acceptance of the detector (cut in theta)
2303     if (fApplyAccCut) { if( partTheta < TMath::Pi()/4.0 || partTheta > 3.0*TMath::Pi()/4.0 ) continue;}
2304
2305     Float_t lambdaTheta = 0.;
2306     Float_t bacTheta    = 0.;
2307     Float_t dghtBarV0Theta = 0.;
2308     Float_t dghtMesV0Theta = 0.;
2309     Float_t bacPt       = 0.;
2310     Float_t dghtBarV0Pt = 0.;
2311     Float_t dghtMesV0Pt = 0.;
2312
2313
2314     if ( fAnalysisType == "ESD" ) {
2315       TParticle* xiMC = lMCstack->Particle( iCurrentLabelStack );
2316       if ( xiMC->GetNDaughters() != 2) continue;
2317       if ( xiMC->GetDaughter(0) < 0 )  continue;
2318       if ( xiMC->GetDaughter(1) < 0 )  continue;
2319
2320       TParticle* lDght0ofXi = lMCstack->Particle(  xiMC->GetDaughter(0) );
2321       TParticle* lDght1ofXi = lMCstack->Particle(  xiMC->GetDaughter(1) );
2322
2323       TParticle* lLambda = 0;
2324       TParticle* lBach   = 0;
2325
2326       // Xi - Case 1
2327       if ( lDght0ofXi->GetPdgCode() == lPdgCodeLambda   &&  // Here !
2328            lDght1ofXi->GetPdgCode() == lPdgCodeBach ){      // Here !
2329
2330         lLambda = lDght0ofXi;
2331         lBach   = lDght1ofXi;
2332       }// end if dghter 0 = Lambda and    dghter 1 = Pi-
2333
2334       // Xi - Case 2
2335         else if ( lDght0ofXi->GetPdgCode() == lPdgCodeBach  &&      // Here !
2336                   lDght1ofXi->GetPdgCode() == lPdgCodeLambda ){     // Here !
2337
2338         lBach   = lDght0ofXi;
2339         lLambda = lDght1ofXi;
2340       }//  end if dghter 0 = Pi-  and   dghter 1 = Lambda
2341
2342         // V0 otherwise - Case 3
2343         else continue;
2344
2345       // Check the emission of particle stays within the acceptance of the detector (cut in pt + theta)
2346       if (fApplyAccCut) {
2347         if ( lLambda->Theta() < TMath::Pi()/4.0  ||    lLambda->Theta() > 3.0*TMath::Pi()/4.0 ) continue;
2348         if( lBach->Theta() < TMath::Pi()/4.0    ||    lBach->Theta() > 3.0*TMath::Pi()/4.0 )   continue;
2349         if( lBach->Pt() < 0.150 ) continue; //FIXME : maybe tuned for Xi but not for K- from Omega ...
2350       }
2351
2352       // -  V0 level ... _____________________________________________________________
2353       TParticle* lDghtBarV0 = 0;
2354       TParticle* lDghtMesV0 = 0;
2355
2356       if( lLambda->GetNDaughters() != 2 )  continue;
2357       if( lLambda->GetDaughter(0) < 0 )    continue;
2358       if( lLambda->GetDaughter(1) < 0 )    continue;
2359
2360
2361       TParticle* lDght0ofLambda = lMCstack->Particle(  lLambda->GetDaughter(0) );
2362       TParticle* lDght1ofLambda = lMCstack->Particle(  lLambda->GetDaughter(1) );
2363
2364       // V0 - Case 1
2365       if ( lDght0ofLambda->GetPdgCode() == lPdgCodeDghtBarV0 &&    // Here !
2366            lDght1ofLambda->GetPdgCode() == lPdgCodeDghtMesV0 ) {    // Here !
2367
2368         lDghtBarV0 = lDght0ofLambda;
2369         lDghtMesV0 = lDght1ofLambda;
2370       }// end if dghter 0 = Proton  and   dghter 1 = Pi-
2371
2372       // V0 - Case 2
2373         else if ( lDght0ofLambda->GetPdgCode() == lPdgCodeDghtMesV0  &&     // Here !
2374                   lDght1ofLambda->GetPdgCode() == lPdgCodeDghtBarV0 ) {      // Here !
2375
2376         lDghtMesV0 = lDght0ofLambda;
2377         lDghtBarV0 = lDght1ofLambda;
2378       }//  end if dghter 0 = Pi-  and   dghter 1 = proton
2379
2380       // V0 otherwise - Case 3
2381         else continue;
2382
2383
2384       // Check the emission of particle stays within the acceptance of the detector
2385       if (fApplyAccCut) {
2386         if( lDghtBarV0->Theta() < TMath::Pi()/4.0  ||  lDghtBarV0->Theta() > 3.0*TMath::Pi()/4.0 ) continue;
2387         if( lDghtMesV0->Theta() < TMath::Pi()/4.0  ||  lDghtMesV0->Theta() > 3.0*TMath::Pi()/4.0 ) continue;
2388
2389         if( lDghtBarV0->Pt() < 0.250 ) continue;
2390         if( lDghtMesV0->Pt() < 0.150 ) continue;
2391       }
2392
2393       lambdaTheta = lLambda->Theta();
2394       bacTheta    = lBach->Theta();
2395       dghtBarV0Theta = lDghtBarV0->Theta();
2396       dghtMesV0Theta = lDghtMesV0->Theta();
2397       bacPt       = lBach->Pt();
2398       dghtBarV0Pt = lDghtBarV0->Pt();
2399       dghtMesV0Pt = lDghtMesV0->Pt();
2400
2401     } else if ( fAnalysisType == "AOD") {
2402
2403       AliAODMCParticle *xiMC = (AliAODMCParticle*) arrayMC->At(iCurrentLabelStack);
2404       if ( xiMC->GetNDaughters() != 2) continue;
2405       if ( xiMC->GetDaughter(0) < 0 )  continue;
2406       if ( xiMC->GetDaughter(1) < 0 )  continue;
2407
2408       AliAODMCParticle* lDght0ofXi = (AliAODMCParticle*) arrayMC->At(  xiMC->GetDaughter(0) );
2409       AliAODMCParticle* lDght1ofXi = (AliAODMCParticle*) arrayMC->At(  xiMC->GetDaughter(1) );
2410
2411       AliAODMCParticle* lLambda = 0;
2412       AliAODMCParticle* lBach   = 0;
2413
2414       // Xi - Case 1
2415       if ( lDght0ofXi->PdgCode() == lPdgCodeLambda   &&  // Here !
2416            lDght1ofXi->PdgCode() == lPdgCodeBach ){      // Here !
2417
2418         lLambda = lDght0ofXi;
2419         lBach   = lDght1ofXi;
2420       }// end if dghter 0 = Lambda and    dghter 1 = Pi-
2421
2422       // Xi - Case 2
2423         else if ( lDght0ofXi->PdgCode() == lPdgCodeBach  &&      // Here !
2424                   lDght1ofXi->PdgCode() == lPdgCodeLambda ){     // Here !
2425
2426         lBach   = lDght0ofXi;
2427         lLambda = lDght1ofXi;
2428       }//  end if dghter 0 = Pi-  and   dghter 1 = Lambda
2429
2430         // V0 otherwise - Case 3
2431         else continue;
2432
2433       // Check the emission of particle stays within the acceptance of the detector (cut in pt + theta)
2434       if (fApplyAccCut) {
2435         if ( lLambda->Theta() < TMath::Pi()/4.0  ||    lLambda->Theta() > 3.0*TMath::Pi()/4.0 ) continue;
2436         if( lBach->Theta() < TMath::Pi()/4.0    ||    lBach->Theta() > 3.0*TMath::Pi()/4.0 )   continue;
2437         if( lBach->Pt() < 0.150 ) continue; //FIXME : maybe tuned for Xi but not for K- from Omega ...
2438       }
2439
2440       // -  V0 level ... _____________________________________________________________
2441       AliAODMCParticle* lDghtBarV0 = 0;
2442       AliAODMCParticle* lDghtMesV0 = 0;
2443
2444       if( lLambda->GetNDaughters() != 2 )  continue;
2445       if( lLambda->GetDaughter(0) < 0 )    continue;
2446       if( lLambda->GetDaughter(1) < 0 )    continue;
2447
2448
2449       AliAODMCParticle* lDght0ofLambda = (AliAODMCParticle*) arrayMC->At(  lLambda->GetDaughter(0) );
2450       AliAODMCParticle* lDght1ofLambda = (AliAODMCParticle*) arrayMC->At(  lLambda->GetDaughter(1) );
2451
2452       // V0 - Case 1
2453       if ( lDght0ofLambda->PdgCode() == lPdgCodeDghtBarV0 &&    // Here !
2454            lDght1ofLambda->PdgCode() == lPdgCodeDghtMesV0 ) {    // Here !
2455
2456         lDghtBarV0 = lDght0ofLambda;
2457         lDghtMesV0 = lDght1ofLambda;
2458       }// end if dghter 0 = Proton  and   dghter 1 = Pi-
2459
2460       // V0 - Case 2
2461         else if ( lDght0ofLambda->PdgCode() == lPdgCodeDghtMesV0  &&     // Here !
2462                   lDght1ofLambda->PdgCode() == lPdgCodeDghtBarV0 ) {      // Here !
2463
2464         lDghtMesV0 = lDght0ofLambda;
2465         lDghtBarV0 = lDght1ofLambda;
2466       }//  end if dghter 0 = Pi-  and   dghter 1 = proton
2467
2468       // V0 otherwise - Case 3
2469         else continue;
2470
2471
2472       // Check the emission of particle stays within the acceptance of the detector
2473       if (fApplyAccCut) {
2474         if( lDghtBarV0->Theta() < TMath::Pi()/4.0  ||  lDghtBarV0->Theta() > 3.0*TMath::Pi()/4.0 ) continue;
2475         if( lDghtMesV0->Theta() < TMath::Pi()/4.0  ||  lDghtMesV0->Theta() > 3.0*TMath::Pi()/4.0 ) continue;
2476
2477         if( lDghtBarV0->Pt() < 0.250 ) continue;
2478         if( lDghtMesV0->Pt() < 0.150 ) continue;
2479       }
2480
2481
2482       lambdaTheta = lLambda->Theta();
2483       bacTheta    = lBach->Theta();
2484       dghtBarV0Theta = lDghtBarV0->Theta();
2485       dghtMesV0Theta = lDghtMesV0->Theta();
2486       bacPt       = lBach->Pt();
2487       dghtBarV0Pt = lDghtBarV0->Pt();
2488       dghtMesV0Pt = lDghtMesV0->Pt();
2489
2490     }
2491
2492     // - Just to know which file is currently open : locate the file containing Xi
2493                 //cout << "Name of the file containing generated Xi :" <<  CurrentFileName()
2494                 //
2495     // - Filling histos for findable cascades... _________________________________________________________________
2496
2497     // - Fill theta histos for Lambda and Bach
2498     lHistThetaLambda        ->Fill( lRadToDeg * lambdaTheta );
2499     lHistThetaBach          ->Fill( lRadToDeg * bacTheta );
2500
2501     // - Fill theta histos for V0 daughters
2502     lHistThetaBarDghter     ->Fill( lRadToDeg * dghtBarV0Theta );
2503     lHistThetaMesDghter     ->Fill( lRadToDeg * dghtMesV0Theta );
2504
2505     // - Fill pt histos.
2506     lHistPtBach             ->Fill( bacPt );
2507     lHistPtBarDghter        ->Fill( dghtBarV0Pt );
2508     lHistPtMesDghter        ->Fill( dghtMesV0Pt );
2509     //  if(iCascType == 1) Printf("Xi- current index = %n ", iCurrentLabelStack);
2510
2511     l2dHistGenPtVsGenYFdbl  ->Fill( partPt, lRapXiMC );
2512
2513     ncascperev++;
2514
2515
2516   }// This is the end of the loop on primaries
2517
2518   if (iCascType == 1) {
2519     fHistnXiMinusPerEv->Fill(ncascperev);
2520     fHistnXiMinusPerEvTot->Fill(ncascperevtot);
2521 //       Printf("N xi-tot = %n N xi-acc = %n\n", ncascperevtot, ncascperev);
2522   }
2523   if (iCascType == 2) {
2524     fHistnXiPlusPerEv->Fill(ncascperev);
2525     fHistnXiPlusPerEvTot->Fill(ncascperevtot);
2526   }
2527   if (iCascType == 3) {
2528     fHistnOmegaMinusPerEv->Fill(ncascperev);
2529     fHistnOmegaMinusPerEvTot->Fill(ncascperevtot);
2530   }
2531   if (iCascType == 4) {
2532     fHistnOmegaPlusPerEv->Fill(ncascperev);
2533     fHistnOmegaPlusPerEvTot->Fill(ncascperevtot);
2534   }
2535
2536
2537 // - Re-initialisation of the local TH1F pointers
2538 lHistEtaGenCasc         = 0x0;
2539 l3dHistGenPtVsGenYvsCentNat = 0x0;
2540 l3dHistGenPtVsGenYvsNtracksNat = 0x0;
2541 l3dHistGenPtVsGenYvsCentInj = 0x0;
2542 l3dHistGenPtVsGenYvsNtracksInj = 0x0;
2543 l3dHistGenPtVsGenctauvsYNat = 0x0;
2544 l3dHistGenPtVsGenctauvsYInj = 0x0;
2545
2546 lHistThetaGenCascNat    = 0x0;
2547 lHistThetaGenCascInj    = 0x0;
2548 l2dHistGenPtVsGenYFdbl  = 0x0;
2549 lHistThetaLambda        = 0x0;
2550 lHistThetaBach          = 0x0;
2551 lHistThetaBarDghter     = 0x0;
2552 lHistThetaMesDghter     = 0x0;
2553 lHistPtBach             = 0x0;
2554 lHistPtBarDghter        = 0x0;
2555 lHistPtMesDghter        = 0x0;
2556
2557 } // end of loop over the different types of cascades (Xi-+, Omega-+)
2558
2559
2560
2561 //__________________________________________________________________________
2562 // Part 3 - Loop over the reconstructed candidates
2563
2564 Int_t nAssoXiMinus = 0;
2565 Int_t nAssoXiPlus = 0;
2566 Int_t nAssoOmegaMinus = 0;
2567 Int_t nAssoOmegaPlus = 0;
2568
2569
2570 Int_t lPosTPCClusters   = 0;
2571 Int_t lNegTPCClusters   = 0;
2572 Int_t lBachTPCClusters  = 0;
2573
2574
2575 // Double_t lChi2Xi         = -1. ;
2576 Double_t lDcaXiDaughters            = -1. ;
2577 Double_t lDcaBachToPrimVertexXi     = -1. ;
2578 Double_t lXiCosineOfPointingAngle   = -1. ;
2579 Double_t lPosXi[3]                  = { -1000.0, -1000.0, -1000.0 };
2580 Double_t lXiRadius                  = -1000. ;
2581
2582 Double_t lInvMassLambdaAsCascDghter = 0.;
2583 Double_t lDcaV0DaughtersXi          = -1.;
2584 // Double_t lV0Chi2Xi               = -1. ;
2585 Double_t lV0CosineOfPointingAngle = -1.;
2586 Double_t lV0toXiCosineOfPointingAngle   = -1.;
2587 Double_t lPosV0Xi[3]                = { -1000. , -1000., -1000. }; // Position of VO coming from cascade
2588 Double_t lV0RadiusXi                = -1000.;
2589 Double_t lDcaV0ToPrimVertexXi       = -1.;
2590 Double_t lDcaPosToPrimVertexXi      = -1.;
2591 Double_t lDcaNegToPrimVertexXi      = -1.;
2592
2593 Double_t lChargeXi                  = -1.;
2594
2595 Double_t lV0mom = -1000.;
2596
2597 Double_t lmcPt = -1.;
2598 Double_t lmcRapCasc = -1.;
2599 Double_t lmcEta = -1000.;
2600 Double_t lmcTransvRadius = -1000.;
2601
2602
2603 Double_t lrecoPt = -100.;
2604 Double_t lrecoTransvRadius = -1000.;
2605
2606 Double_t lDeltaPhiMcReco = -1.;
2607
2608 Double_t lmcPtPosV0Dghter = -100.;
2609 Double_t lmcPtNegV0Dghter = -100.;
2610 Double_t lrecoP = -100.;
2611
2612 Double_t   lmcPtBach = -100.;
2613
2614 Double_t cascadeMass = 0.;
2615
2616 for (Int_t iXi = 0; iXi < ncascades; iXi++) {// This is the begining of the Cascade loop
2617
2618
2619       Bool_t   lIsPosInXiProton      = kFALSE;
2620       Bool_t   lIsPosInXiPion        = kFALSE;
2621       Bool_t   lIsPosInOmegaProton   = kFALSE;
2622       Bool_t   lIsPosInOmegaPion     = kFALSE;
2623
2624       Bool_t   lIsNegInXiProton      = kFALSE;
2625       Bool_t   lIsNegInXiPion        = kFALSE;
2626       Bool_t   lIsNegInOmegaProton   = kFALSE;
2627       Bool_t   lIsNegInOmegaPion     = kFALSE;
2628
2629       Bool_t   lIsBachelorKaon       = kFALSE;
2630       Bool_t   lIsBachelorPion       = kFALSE;
2631
2632       Bool_t   lIsBachelorKaonForTPC = kFALSE;
2633       Bool_t   lIsBachelorPionForTPC = kFALSE;
2634       Bool_t   lIsNegPionForTPC      = kFALSE;
2635       Bool_t   lIsPosPionForTPC      = kFALSE;
2636       Bool_t   lIsNegProtonForTPC    = kFALSE;
2637       Bool_t   lIsPosProtonForTPC    = kFALSE;
2638
2639       // Combined PID
2640       // Reasonable guess for the priors for the cascade track sample (e-, mu, pi, K, p)
2641       Double_t lPriorsGuessXi[5]    = {0, 0, 2, 0, 1};
2642       Double_t lPriorsGuessOmega[5] = {0, 0, 1, 1, 1};
2643
2644       Double_t ppionBach = 0.0, pkaonBach = 0.0;
2645
2646       Bool_t   lIsBachelorMCPiMinus  = kFALSE;
2647       Bool_t   lIsBachelorMCPiPlus   = kFALSE;
2648       Bool_t   lIsBachelorMCKMinus   = kFALSE;
2649       Bool_t   lIsBachelorMCKPlus    = kFALSE;
2650
2651       Double_t lInvMassXiMinus    = 0.;
2652       Double_t lInvMassXiPlus     = 0.;
2653       Double_t lInvMassOmegaMinus = 0.;
2654       Double_t lInvMassOmegaPlus  = 0.;
2655
2656       Bool_t lAssoXiMinus    = kFALSE;
2657       Bool_t lAssoXiPlus     = kFALSE;
2658       Bool_t lAssoOmegaMinus = kFALSE;
2659       Bool_t lAssoOmegaPlus  = kFALSE;
2660
2661       Bool_t kIsNaturalPart = kFALSE;
2662
2663       if ( fAnalysisType == "ESD" ) {
2664         AliESDcascade *xiESD = lESDevent->GetCascade(iXi);
2665         if (!xiESD) continue;
2666
2667         // - Step II.1 : Connection to daughter tracks of the current cascade
2668         //-------------
2669
2670                 UInt_t lIdxPosXi        = (UInt_t) TMath::Abs( xiESD->GetPindex() );
2671                 UInt_t lIdxNegXi        = (UInt_t) TMath::Abs( xiESD->GetNindex() );
2672                 UInt_t lBachIdx         = (UInt_t) TMath::Abs( xiESD->GetBindex() );
2673                 // abs value not needed ; the index should always be positive (!= label ...)
2674
2675
2676         // FIXME : rejection of a double use of a daughter track (nothing but just a crosscheck of what is done in the cascade vertexer)
2677         if(lBachIdx == lIdxNegXi) {
2678                 AliWarning("Pb / Idx(Bach. track) = Idx(Neg. track) ... continue!"); continue;
2679         }
2680         if(lBachIdx == lIdxPosXi) {
2681                 AliWarning("Pb / Idx(Bach. track) = Idx(Pos. track) ... continue!"); continue;
2682         }
2683
2684         AliESDtrack *pTrackXi           = lESDevent->GetTrack( lIdxPosXi );
2685         AliESDtrack *nTrackXi           = lESDevent->GetTrack( lIdxNegXi );
2686         AliESDtrack *bachTrackXi        = lESDevent->GetTrack( lBachIdx  );
2687         if (!pTrackXi || !nTrackXi || !bachTrackXi ) {
2688                 Printf("ERROR: Could not retrieve one of the 3 daughter tracks of the cascade ...");
2689                 continue;
2690         }
2691
2692         lPosTPCClusters   = pTrackXi->GetTPCNcls();
2693         lNegTPCClusters   = nTrackXi->GetTPCNcls();
2694         lBachTPCClusters  = bachTrackXi->GetTPCNcls();
2695
2696                 // FIXME : rejection of a poor quality tracks
2697         if(fkQualityCutTPCrefit){
2698                 // 1 - Poor quality related to TPCrefit
2699                 ULong_t pStatus    = pTrackXi->GetStatus();
2700                 ULong_t nStatus    = nTrackXi->GetStatus();
2701                 ULong_t bachStatus = bachTrackXi->GetStatus();
2702                 if ((pStatus&AliESDtrack::kTPCrefit)    == 0) { AliWarning("Pb / V0 Pos. track has no TPCrefit ... continue!"); continue; }
2703                 if ((nStatus&AliESDtrack::kTPCrefit)    == 0) { AliWarning("Pb / V0 Neg. track has no TPCrefit ... continue!"); continue; }
2704                 if ((bachStatus&AliESDtrack::kTPCrefit) == 0) { AliWarning("Pb / Bach.   track has no TPCrefit ... continue!"); continue; }
2705         }
2706         if(fkQualityCutnTPCcls){
2707                 // 2 - Poor quality related to TPC clusters
2708                 if(lPosTPCClusters  < fMinnTPCcls) { AliWarning("Pb / V0 Pos. track has less than n TPC clusters ... continue!"); continue; }
2709                 if(lNegTPCClusters  < fMinnTPCcls) { AliWarning("Pb / V0 Neg. track has less than n TPC clusters ... continue!"); continue; }
2710                 if(lBachTPCClusters < fMinnTPCcls) { AliWarning("Pb / Bach.   track has less than n TPC clusters ... continue!"); continue; }
2711         }
2712
2713         // - Step II.2 : Info over reconstructed cascades
2714         //-------------
2715
2716
2717         Double_t lV0quality = 0.;
2718
2719         if( bachTrackXi->Charge() < 0 ) {
2720                 lV0quality = 0.;
2721                 xiESD->ChangeMassHypothesis(lV0quality , 3312);
2722                         // Calculate the effective mass of the Xi- candidate.
2723                         // pdg code 3312 = Xi-
2724                 lInvMassXiMinus = xiESD->GetEffMassXi();
2725
2726                 lV0quality = 0.;
2727                 xiESD->ChangeMassHypothesis(lV0quality , 3334);
2728                         // Calculate the effective mass of the Xi- candidate.
2729                         // pdg code 3334 = Omega-
2730                 lInvMassOmegaMinus = xiESD->GetEffMassXi();
2731
2732                 lV0quality = 0.;
2733                 xiESD->ChangeMassHypothesis(lV0quality , 3312);         // Back to default hyp.
2734
2735         }
2736
2737         if( bachTrackXi->Charge() >  0 ){
2738                 lV0quality = 0.;
2739                 xiESD->ChangeMassHypothesis(lV0quality , -3312);
2740                         // Calculate the effective mass of the Xi+ candidate.
2741                         // pdg code -3312 = Xi+
2742                 lInvMassXiPlus = xiESD->GetEffMassXi();
2743
2744                 lV0quality = 0.;
2745                 xiESD->ChangeMassHypothesis(lV0quality , -3334);
2746                         // Calculate the effective mass of the Xi+ candidate.
2747                         // pdg code -3334  = Omega+
2748                 lInvMassOmegaPlus = xiESD->GetEffMassXi();
2749
2750                 lV0quality = 0.;
2751                 xiESD->ChangeMassHypothesis(lV0quality , -3312);        // Back to "default" hyp.
2752         }
2753
2754
2755         //lChi2Xi                       = xiESD->GetChi2Xi();
2756         lDcaXiDaughters                 = xiESD->GetDcaXiDaughters();
2757         lDcaBachToPrimVertexXi          = TMath::Abs( bachTrackXi->GetD( lBestPrimaryVtxPos[0],
2758                                                                          lBestPrimaryVtxPos[1],
2759                                                                          lMagneticField  ) );
2760                                         // NOTE : AliExternalTrackParam::GetD returns an algebraic value
2761         lXiCosineOfPointingAngle        = xiESD->GetCascadeCosineOfPointingAngle( lBestPrimaryVtxPos[0],
2762                                                                                   lBestPrimaryVtxPos[1],
2763                                                                                   lBestPrimaryVtxPos[2] );
2764                                         // Take care : the best available vertex should be used (like in AliCascadeVertexer)
2765         xiESD->GetXYZcascade( lPosXi[0],  lPosXi[1], lPosXi[2] );
2766         lInvMassLambdaAsCascDghter      = xiESD->GetEffMass();
2767                                         // This value shouldn't change, whatever the working hyp. is : Xi-, Xi+, Omega-, Omega+
2768         lDcaV0DaughtersXi               = xiESD->GetDcaV0Daughters();
2769         // lV0Chi2Xi                    = xiESD->GetChi2V0();
2770         lV0toXiCosineOfPointingAngle    = xiESD->GetV0CosineOfPointingAngle( lPosXi[0],
2771                                                                              lPosXi[1],
2772                                                                              lPosXi[2] );
2773
2774         lV0CosineOfPointingAngle      = xiESD->GetV0CosineOfPointingAngle( lBestPrimaryVtxPos[0],
2775                                                                            lBestPrimaryVtxPos[1],
2776                                                                            lBestPrimaryVtxPos[2]);
2777
2778         //if (lV0toXiCosineOfPointingAngle==1.) cout << "Cosine V0 PA wrt Xi pos ==1!" <<endl;
2779         //else if (lV0toXiCosineOfPointingAngle>1.) cout <<"Cosine V0 PA wrt Xi pos >1!"<<endl;
2780         xiESD->GetXYZ( lPosV0Xi[0],  lPosV0Xi[1], lPosV0Xi[2] );
2781
2782         lDcaV0ToPrimVertexXi            = xiESD->GetD( lBestPrimaryVtxPos[0],
2783                                                 lBestPrimaryVtxPos[1],
2784                                                 lBestPrimaryVtxPos[2] );
2785
2786         lDcaPosToPrimVertexXi           = TMath::Abs( pTrackXi   ->GetD( lBestPrimaryVtxPos[0],
2787                                                                          lBestPrimaryVtxPos[1],
2788                                                                          lMagneticField  )     );
2789
2790         lDcaNegToPrimVertexXi           = TMath::Abs( nTrackXi   ->GetD( lBestPrimaryVtxPos[0],
2791                                                                          lBestPrimaryVtxPos[1],
2792                                                                          lMagneticField  )     );
2793
2794         lChargeXi = xiESD->Charge();
2795
2796         // - Step II.3 : PID info
2797
2798         //-------------
2799
2800
2801         // 3.1 - PID Information
2802
2803
2804         // Combined VO-positive-daughter PID
2805         AliPID pPidXi;          pPidXi.SetPriors(    lPriorsGuessXi    );
2806         AliPID pPidOmega;       pPidOmega.SetPriors( lPriorsGuessOmega );
2807
2808         if( pTrackXi->IsOn(AliESDtrack::kESDpid) ){  // Combined PID exists
2809                 Double_t r[10] = {0.}; pTrackXi->GetESDpid(r);
2810                 pPidXi.SetProbabilities(r);
2811                 pPidOmega.SetProbabilities(r);
2812
2813                 // Check if the V0 positive track is a proton (case for Xi-)
2814                 Double_t pproton = pPidXi.GetProbability(AliPID::kProton);
2815                 if (pproton > pPidXi.GetProbability(AliPID::kElectron) &&
2816                     pproton > pPidXi.GetProbability(AliPID::kMuon)     &&
2817                     pproton > pPidXi.GetProbability(AliPID::kPion)     &&
2818                     pproton > pPidXi.GetProbability(AliPID::kKaon)     )     lIsPosInXiProton = kTRUE;
2819
2820                 // Check if the V0 positive track is a pi+ (case for Xi+)
2821                 Double_t ppion = pPidXi.GetProbability(AliPID::kPion);
2822                 if (ppion > pPidXi.GetProbability(AliPID::kElectron) &&
2823                     ppion > pPidXi.GetProbability(AliPID::kMuon)     &&
2824                     ppion > pPidXi.GetProbability(AliPID::kKaon)     &&
2825                     ppion > pPidXi.GetProbability(AliPID::kProton)   )     lIsPosInXiPion = kTRUE;
2826
2827
2828                 // Check if the V0 positive track is a proton (case for Omega-)
2829                 pproton = 0.;
2830                     pproton = pPidOmega.GetProbability(AliPID::kProton);
2831                 if (pproton > pPidOmega.GetProbability(AliPID::kElectron) &&
2832                     pproton > pPidOmega.GetProbability(AliPID::kMuon)     &&
2833                     pproton > pPidOmega.GetProbability(AliPID::kPion)     &&
2834                     pproton > pPidOmega.GetProbability(AliPID::kKaon)     )  lIsPosInOmegaProton = kTRUE;
2835
2836                 // Check if the V0 positive track is a pi+ (case for Omega+)
2837                 ppion = 0.;
2838                     ppion = pPidOmega.GetProbability(AliPID::kPion);
2839                 if (ppion > pPidOmega.GetProbability(AliPID::kElectron) &&
2840                     ppion > pPidOmega.GetProbability(AliPID::kMuon)     &&
2841                     ppion > pPidOmega.GetProbability(AliPID::kKaon)     &&
2842                     ppion > pPidOmega.GetProbability(AliPID::kProton)   )    lIsPosInOmegaPion = kTRUE;
2843
2844         }// end if V0 positive track with existing combined PID
2845
2846
2847         // Combined VO-negative-daughter PID
2848         AliPID nPidXi;          nPidXi.SetPriors(    lPriorsGuessXi    );
2849         AliPID nPidOmega;       nPidOmega.SetPriors( lPriorsGuessOmega );
2850
2851         if( nTrackXi->IsOn(AliESDtrack::kESDpid) ) {  // Combined PID exists
2852                 Double_t r[10] = {0.}; nTrackXi->GetESDpid(r);
2853                 nPidXi.SetProbabilities(r);
2854                 nPidOmega.SetProbabilities(r);
2855
2856                 // Check if the V0 negative track is a pi- (case for Xi-)
2857                 Double_t ppion = nPidXi.GetProbability(AliPID::kPion);
2858                 if (ppion > nPidXi.GetProbability(AliPID::kElectron) &&
2859                     ppion > nPidXi.GetProbability(AliPID::kMuon)     &&
2860                     ppion > nPidXi.GetProbability(AliPID::kKaon)     &&
2861                     ppion > nPidXi.GetProbability(AliPID::kProton)   )     lIsNegInXiPion = kTRUE;
2862
2863                 // Check if the V0 negative track is an anti-proton (case for Xi+)
2864                 Double_t pproton = nPidXi.GetProbability(AliPID::kProton);
2865                 if (pproton > nPidXi.GetProbability(AliPID::kElectron) &&
2866                     pproton > nPidXi.GetProbability(AliPID::kMuon)     &&
2867                     pproton > nPidXi.GetProbability(AliPID::kPion)     &&
2868                     pproton > nPidXi.GetProbability(AliPID::kKaon)     )     lIsNegInXiProton = kTRUE;
2869
2870                 // Check if the V0 negative track is a pi- (case for Omega-)
2871                 ppion = 0.;
2872                     ppion = nPidOmega.GetProbability(AliPID::kPion);
2873                 if (ppion > nPidOmega.GetProbability(AliPID::kElectron) &&
2874                     ppion > nPidOmega.GetProbability(AliPID::kMuon)     &&
2875                     ppion > nPidOmega.GetProbability(AliPID::kKaon)     &&
2876                     ppion > nPidOmega.GetProbability(AliPID::kProton)   )    lIsNegInOmegaPion = kTRUE;
2877
2878                 // Check if the V0 negative track is an anti-proton (case for Omega+)
2879                 pproton = 0.;
2880                     pproton = nPidOmega.GetProbability(AliPID::kProton);
2881                 if (pproton > nPidOmega.GetProbability(AliPID::kElectron) &&
2882                     pproton > nPidOmega.GetProbability(AliPID::kMuon)     &&
2883                     pproton > nPidOmega.GetProbability(AliPID::kPion)     &&
2884                     pproton > nPidOmega.GetProbability(AliPID::kKaon)     )  lIsNegInOmegaProton = kTRUE;
2885
2886         }// end if V0 negative track with existing combined PID
2887
2888
2889         // Combined bachelor PID
2890         AliPID bachPidXi;       bachPidXi.SetPriors(    lPriorsGuessXi    );
2891         AliPID bachPidOmega;    bachPidOmega.SetPriors( lPriorsGuessOmega );
2892
2893
2894         if ( bachTrackXi->IsOn(AliESDtrack::kESDpid) ) {  // Combined PID exists
2895                 Double_t r[10] = {0.}; bachTrackXi->GetESDpid(r);
2896                 bachPidXi.SetProbabilities(r);
2897                 bachPidOmega.SetProbabilities(r);
2898                 // Check if the bachelor track is a pion
2899                     ppionBach = bachPidXi.GetProbability(AliPID::kPion);
2900                 if (ppionBach > bachPidXi.GetProbability(AliPID::kElectron) &&
2901                     ppionBach > bachPidXi.GetProbability(AliPID::kMuon)     &&
2902                     ppionBach > bachPidXi.GetProbability(AliPID::kKaon)     &&
2903                     ppionBach > bachPidXi.GetProbability(AliPID::kProton)   )     lIsBachelorPion = kTRUE;
2904                 // Check if the bachelor track is a kaon
2905                     pkaonBach = bachPidOmega.GetProbability(AliPID::kKaon);
2906                 if (pkaonBach > bachPidOmega.GetProbability(AliPID::kElectron) &&
2907                     pkaonBach > bachPidOmega.GetProbability(AliPID::kMuon)     &&
2908                     pkaonBach > bachPidOmega.GetProbability(AliPID::kPion)     &&
2909                     pkaonBach > bachPidOmega.GetProbability(AliPID::kProton)   )  lIsBachelorKaon = kTRUE;
2910         }// end if bachelor track with existing combined PID
2911
2912
2913         // 3.1.B - TPC PID : 4-sigma bands on Bethe-Bloch curve
2914
2915         // Bachelor
2916         if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( bachTrackXi,AliPID::kKaon)) < 4) lIsBachelorKaonForTPC = kTRUE;
2917         if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( bachTrackXi,AliPID::kPion)) < 4) lIsBachelorPionForTPC = kTRUE;
2918
2919         // Negative V0 daughter
2920         if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( nTrackXi,AliPID::kPion   )) < 4) lIsNegPionForTPC   = kTRUE;
2921         if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( nTrackXi,AliPID::kProton )) < 4) lIsNegProtonForTPC = kTRUE;
2922
2923         // Positive V0 daughter
2924         if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( pTrackXi,AliPID::kPion   )) < 4) lIsPosPionForTPC   = kTRUE;
2925         if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( pTrackXi,AliPID::kProton )) < 4) lIsPosProtonForTPC = kTRUE;
2926
2927 /*
2928         const AliExternalTrackParam *pInnerWallTrackXi    = pTrackXi    ->GetInnerParam(); // Do not use GetTPCInnerWall
2929         const AliExternalTrackParam *nInnerWallTrackXi    = nTrackXi    ->GetInnerParam();
2930         const AliExternalTrackParam *bachInnerWallTrackXi = bachTrackXi ->GetInnerParam();
2931         if(pInnerWallTrackXi && nInnerWallTrackXi && bachInnerWallTrackXi ){
2932
2933                 Double_t pMomInnerWall    = pInnerWallTrackXi   ->GetP();
2934                 Double_t nMomInnerWall    = nInnerWallTrackXi   ->GetP();
2935                 Double_t bachMomInnerWall = bachInnerWallTrackXi->GetP();
2936
2937                 // Bachelor
2938                 if (TMath::Abs(fESDpid->NumberOfSigmasTPC( bachTrackXi,AliPID::kPion)) < 3)                              lIsBachelorPionForTPC = kTRUE;
2939                 if (bachMomInnerWall < 0.350  && TMath::Abs(fESDpid->NumberOfSigmasTPC( bachTrackXi,AliPID::kKaon)) < 5) lIsBachelorKaonForTPC = kTRUE;
2940                 if (bachMomInnerWall > 0.350  && TMath::Abs(fESDpid->NumberOfSigmasTPC( bachTrackXi,AliPID::kKaon)) < 3) lIsBachelorKaonForTPC = kTRUE;
2941
2942                 // Negative V0 daughter
2943                 if (TMath::Abs(fESDpid->NumberOfSigmasTPC( nTrackXi,AliPID::kPion   )) < 3  )                           lIsNegPionForTPC   = kTRUE;
2944                 if (nMomInnerWall < 0.6  && TMath::Abs(fESDpid->NumberOfSigmasTPC( nTrackXi,AliPID::kProton ) ) < 5 )   lIsNegProtonForTPC = kTRUE;
2945                 if (nMomInnerWall > 0.6  && TMath::Abs(fESDpid->NumberOfSigmasTPC( nTrackXi,AliPID::kProton ) ) < 3 )   lIsNegProtonForTPC = kTRUE;
2946
2947                 // Positive V0 daughter
2948                 if (TMath::Abs(fESDpid->NumberOfSigmasTPC( pTrackXi,AliPID::kPion   )) < 3 )                            lIsPosPionForTPC   = kTRUE;
2949                 if (pMomInnerWall < 0.6  && TMath::Abs(fESDpid->NumberOfSigmasTPC( pTrackXi,AliPID::kProton )) < 5)     lIsPosProtonForTPC = kTRUE;
2950                 if (pMomInnerWall > 0.6  && TMath::Abs(fESDpid->NumberOfSigmasTPC( pTrackXi,AliPID::kProton )) < 3)     lIsPosProtonForTPC = kTRUE;
2951
2952         }
2953 */
2954
2955         // 3.2 - PID proba Vs Pt(Bach)
2956         Int_t      lblBachForPID  = (Int_t) TMath::Abs( bachTrackXi->GetLabel() );
2957         TParticle* mcBachForPID   = lMCstack->Particle( lblBachForPID );
2958         lmcPtBach      = mcBachForPID->Pt();
2959
2960         // 3.3 - MC perfect PID
2961
2962         if( mcBachForPID->GetPdgCode() == -211) lIsBachelorMCPiMinus = kTRUE;
2963         if( mcBachForPID->GetPdgCode() ==  211) lIsBachelorMCPiPlus  = kTRUE;
2964         if( mcBachForPID->GetPdgCode() == -321) lIsBachelorMCKMinus  = kTRUE;
2965         if( mcBachForPID->GetPdgCode() ==  321) lIsBachelorMCKPlus   = kTRUE;
2966
2967
2968
2969         // - Step II.4 : MC association (care : lots of "continue;" below this line)
2970         //-------------
2971
2972         if(fDebug > 5)
2973                 cout    << "MC EventNumber : " << lMCevent->Header()->GetEvent()
2974                         << " / MC event Number in Run : " << lMCevent->Header()->GetEventNrInRun() << endl;
2975
2976         // - Step 4.1 : level of the V0 daughters
2977
2978         Int_t lblPosV0Dghter = (Int_t) TMath::Abs( pTrackXi->GetLabel() );
2979                 // Abs value = needed ! question of quality track association ... (negative label when at least one cluster in the track is from a different particle)
2980         Int_t lblNegV0Dghter = (Int_t) TMath::Abs( nTrackXi->GetLabel() );
2981         TParticle* mcPosV0Dghter = lMCstack->Particle( lblPosV0Dghter );
2982         TParticle* mcNegV0Dghter = lMCstack->Particle( lblNegV0Dghter );
2983
2984
2985         // - Step 4.2 : level of the Xi daughters
2986
2987         Int_t lblMotherPosV0Dghter = mcPosV0Dghter->GetFirstMother() ;
2988         Int_t lblMotherNegV0Dghter = mcNegV0Dghter->GetFirstMother();
2989
2990                 if( lblMotherPosV0Dghter != lblMotherNegV0Dghter) continue; // same mother
2991                 if( lblMotherPosV0Dghter < 0 ) continue; // this particle is primary, no mother
2992                 if( lblMotherNegV0Dghter < 0 ) continue;
2993
2994
2995                 // mothers = Lambda candidate ... a priori
2996
2997         TParticle* mcMotherPosV0Dghter = lMCstack->Particle( lblMotherPosV0Dghter );
2998         TParticle* mcMotherNegV0Dghter = lMCstack->Particle( lblMotherNegV0Dghter );  // MN: redundant?? already checked that labels are the same...-->same part from stack
2999
3000         Int_t      lblBach  = (Int_t) TMath::Abs( bachTrackXi->GetLabel() );
3001         TParticle* mcBach   = lMCstack->Particle( lblBach );
3002
3003
3004         // - Step 4.3 : level of Xi candidate
3005
3006         Int_t lblGdMotherPosV0Dghter =   mcMotherPosV0Dghter->GetFirstMother() ;
3007         Int_t lblGdMotherNegV0Dghter =   mcMotherNegV0Dghter->GetFirstMother() ;
3008
3009                 if( lblGdMotherPosV0Dghter != lblGdMotherNegV0Dghter ) continue;
3010                 if( lblGdMotherPosV0Dghter < 0 ) continue; // primary lambda ...
3011                 if( lblGdMotherNegV0Dghter < 0 ) continue; // primary lambda ...
3012
3013                 // Gd mothers = Xi candidate ... a priori
3014
3015         TParticle* mcGdMotherPosV0Dghter = lMCstack->Particle( lblGdMotherPosV0Dghter );
3016         TParticle* mcGdMotherNegV0Dghter = lMCstack->Particle( lblGdMotherNegV0Dghter );
3017
3018         Int_t lblMotherBach = (Int_t) TMath::Abs( mcBach->GetFirstMother()  );
3019
3020         if( lblMotherBach != lblGdMotherPosV0Dghter ) continue; //same mother for bach and V0 daughters
3021
3022         TParticle* mcMotherBach = lMCstack->Particle( lblMotherBach );
3023
3024         // Check if cascade is primary
3025
3026         if (!(lMCstack->IsPhysicalPrimary(lblMotherBach))) continue;
3027
3028         kIsNaturalPart = lMCevent->IsFromBGEvent(lblMotherBach); // if kTRUE the particle is a natural one
3029         if (!kIsNaturalPart) {
3030           if (!(mcMotherBach->GetFirstMother()<0)) kIsNaturalPart = kTRUE;
3031         }
3032
3033         // - Step 4.4 : Manage boolean for association
3034
3035         if( mcMotherBach                ->GetPdgCode() ==   3312 &&
3036             mcGdMotherPosV0Dghter       ->GetPdgCode() ==   3312 &&
3037             mcGdMotherNegV0Dghter       ->GetPdgCode() ==   3312)    {  lAssoXiMinus = kTRUE;
3038                                                                         cascadeMass = 1.321;
3039                                                                       //  Printf("Xi- asso current index = %n ", lblGdMotherPosV0Dghter);
3040                                                                         nAssoXiMinus++; }
3041
3042         else if( mcMotherBach           ->GetPdgCode() ==  -3312 &&
3043             mcGdMotherPosV0Dghter       ->GetPdgCode() ==  -3312 &&
3044             mcGdMotherNegV0Dghter       ->GetPdgCode() ==  -3312)    {  lAssoXiPlus = kTRUE;
3045                                                                         cascadeMass = 1.321;
3046                                                                         nAssoXiPlus++; }
3047
3048         else if( mcMotherBach           ->GetPdgCode() ==   3334 &&
3049             mcGdMotherPosV0Dghter       ->GetPdgCode() ==   3334 &&
3050             mcGdMotherNegV0Dghter       ->GetPdgCode() ==   3334)    {  lAssoOmegaMinus = kTRUE;
3051                                                                         cascadeMass = 1.672;
3052                                                                         nAssoOmegaMinus++; }
3053
3054         else if( mcMotherBach           ->GetPdgCode() ==  -3334 &&
3055             mcGdMotherPosV0Dghter       ->GetPdgCode() ==  -3334 &&
3056             mcGdMotherNegV0Dghter       ->GetPdgCode() ==  -3334)    {  lAssoOmegaPlus = kTRUE;
3057                                                                         cascadeMass = 1.672;
3058                                                                         nAssoOmegaPlus++; }
3059
3060
3061
3062         // If a proper association  exists ...
3063
3064         if(fDebug > 4){
3065                 cout << "XiMinus    = " << lAssoXiMinus    << endl;
3066                 cout << "XiPlus     = " << lAssoXiPlus     << endl;
3067                 cout << "OmegaMinus = " << lAssoOmegaMinus << endl;
3068                 cout << "OmegaPlus  = " << lAssoOmegaPlus  << endl
3069                      << "----"          << endl;
3070         }
3071
3072
3073         if(fDebug > 5){
3074                 cout << endl;
3075                 cout << "- V0 daughters - " << endl;
3076                 cout << "     + V0 Pos. / Label : " << lblPosV0Dghter
3077                 << " - Pdg Code : " << mcPosV0Dghter->GetTitle() << endl;
3078                 cout << "     - V0 Neg. / Label : " << lblNegV0Dghter
3079                 << " - Pdg Code : " << mcNegV0Dghter->GetTitle() << endl;
3080
3081                 cout << "- Xi daughters - " << endl;
3082                 cout << "     + V0 Pos. mother / Label : " << lblMotherPosV0Dghter
3083                 << " - Pdg Code : " << mcMotherPosV0Dghter->GetTitle() << endl;
3084                 cout << "     - V0 Neg. mother / Label : " << lblMotherNegV0Dghter
3085                 << " - Pdg Code : " << mcMotherNegV0Dghter->GetTitle() << endl;
3086
3087                 cout << "     --  Bach. / Label :" << lblBach
3088                 << " -  Pdg Code : " << mcBach->GetTitle() << endl;
3089
3090                 cout << "- Xi candidate -" << endl;
3091                 cout << "    +  V0 Pos. Gd Mother / Label : " << lblGdMotherPosV0Dghter
3092                 << " - Pdg Code : " << mcGdMotherPosV0Dghter->GetTitle() << endl;
3093                 cout << "    -  V0 Neg. Gd Mother / Label : "  << lblGdMotherNegV0Dghter
3094                 << " - Pdg Code : "<< mcGdMotherNegV0Dghter->GetTitle() << endl;
3095
3096                 cout << "    --  Mother Bach. / Label : " << lblMotherBach
3097                 << " - Pdg Code    : " << mcMotherBach->GetTitle() << endl;
3098                 cout << endl;
3099         }
3100
3101
3102         lmcPt             = mcMotherBach->Pt();
3103         lmcRapCasc        = 0.5*TMath::Log( (mcMotherBach->Energy() + mcMotherBach->Pz()) /
3104                                                      (mcMotherBach->Energy() - mcMotherBach->Pz() +1.e-13) );
3105         lmcEta            = mcMotherBach->Eta();
3106         lmcTransvRadius   = mcBach->R(); // to get the decay point of Xi, = the production vertex of Bachelor ...
3107
3108         TVector3 lmcTVect3Mom( mcMotherBach->Px(), mcMotherBach->Py(), mcMotherBach->Pz() );
3109
3110         lrecoPt           = xiESD->Pt();
3111         lrecoTransvRadius = TMath::Sqrt( xiESD->Xv() * xiESD->Xv() + xiESD->Yv() * xiESD->Yv() );
3112
3113         TVector3 lrecoTVect3Mom( xiESD->Px(), xiESD->Py(), xiESD->Pz() );
3114         lDeltaPhiMcReco   = lmcTVect3Mom.DeltaPhi( lrecoTVect3Mom ) * 180.0/TMath::Pi();
3115
3116         lmcPtPosV0Dghter = mcPosV0Dghter->Pt() ;
3117         lmcPtNegV0Dghter = mcNegV0Dghter->Pt();
3118         lrecoP           = xiESD->P();
3119
3120         Double_t nV0mom[3] = {0. ,0. ,0. };
3121         Double_t pV0mom[3] = {0. ,0. ,0. };
3122
3123         xiESD->GetNPxPyPz(nV0mom[0],nV0mom[1],nV0mom[2]);
3124         xiESD->GetPPxPyPz(pV0mom[0],pV0mom[1],pV0mom[2]);
3125
3126         lV0mom = TMath::Sqrt(TMath::Power(nV0mom[0]+pV0mom[0],2)+TMath::Power(nV0mom[1]+pV0mom[1],2)+TMath::Power(nV0mom[2]+pV0mom[2],2));
3127
3128       } else if ( fAnalysisType == "AOD" ) {
3129
3130         const AliAODcascade *xiAOD = lAODevent->GetCascade(iXi);
3131
3132         if (!xiAOD) continue;
3133
3134         // - Step II.1 : Connection to daughter tracks of the current cascade
3135         //-------------
3136
3137         AliAODTrack *pTrackXi    = dynamic_cast<AliAODTrack*>( xiAOD->GetDaughter(0) );
3138         AliAODTrack *nTrackXi    = dynamic_cast<AliAODTrack*>( xiAOD->GetDaughter(1) );
3139         AliAODTrack *bachTrackXi = dynamic_cast<AliAODTrack*>( xiAOD->GetDecayVertexXi()->GetDaughter(0) );
3140         if (!pTrackXi || !nTrackXi || !bachTrackXi ) {
3141           AliWarning("ERROR: Could not retrieve one of the 3 AOD daughter tracks of the cascade ...");
3142           continue;
3143         }
3144
3145         UInt_t lIdxPosXi  = (UInt_t) TMath::Abs( pTrackXi->GetID() );
3146         UInt_t lIdxNegXi  = (UInt_t) TMath::Abs( nTrackXi->GetID() );
3147         UInt_t lBachIdx   = (UInt_t) TMath::Abs( bachTrackXi->GetID() );
3148
3149                 // abs value not needed ; the index should always be positive (!= label ...)
3150
3151
3152         // FIXME : rejection of a double use of a daughter track (nothing but just a crosscheck of what is done in the cascade vertexer)
3153         if(lBachIdx == lIdxNegXi) {
3154                 AliWarning("Pb / Idx(Bach. track) = Idx(Neg. track) ... continue!"); continue;
3155         }
3156         if(lBachIdx == lIdxPosXi) {
3157                 AliWarning("Pb / Idx(Bach. track) = Idx(Pos. track) ... continue!"); continue;
3158         }
3159
3160         lPosTPCClusters   = pTrackXi->GetTPCNcls();
3161         lNegTPCClusters   = nTrackXi->GetTPCNcls();
3162         lBachTPCClusters  = bachTrackXi->GetTPCNcls();
3163
3164
3165         // FIXME : rejection of a poor quality tracks
3166         if (fkQualityCutTPCrefit) {
3167                 // 1 - Poor quality related to TPCrefit
3168           if (!(pTrackXi->IsOn(AliAODTrack::kTPCrefit))) { AliWarning("Pb / V0 Pos. track has no TPCrefit ... continue!"); continue; }
3169           if (!(nTrackXi->IsOn(AliAODTrack::kTPCrefit))) { AliWarning("Pb / V0 Neg. track has no TPCrefit ... continue!"); continue; }
3170           if (!(bachTrackXi->IsOn(AliAODTrack::kTPCrefit))) { AliWarning("Pb / Bach.   track has no TPCrefit ... continue!"); continue; }
3171
3172         }
3173         if (fkQualityCutnTPCcls) {
3174                 // 2 - Poor quality related to TPC clusters
3175                 if(lPosTPCClusters  < fMinnTPCcls) { AliWarning("Pb / V0 Pos. track has less than n TPC clusters ... continue!"); continue; }
3176                 if(lNegTPCClusters  < fMinnTPCcls) { AliWarning("Pb / V0 Neg. track has less than n TPC clusters ... continue!"); continue; }
3177                 if(lBachTPCClusters < fMinnTPCcls) { AliWarning("Pb / Bach.   track has less than n TPC clusters ... continue!"); continue; }
3178         }
3179
3180         // - Step II.2 : Info over reconstructed cascades
3181         //-------------
3182
3183
3184         if( bachTrackXi->Charge() < 0 ) {
3185                 lInvMassXiMinus = xiAOD->MassXi();
3186
3187                 lInvMassOmegaMinus = xiAOD->MassOmega();
3188         }
3189
3190         if( bachTrackXi->Charge() >  0 ){
3191                 lInvMassXiPlus = xiAOD->MassXi();
3192
3193                 lInvMassOmegaPlus = xiAOD->MassOmega();
3194
3195         }
3196
3197
3198         //lChi2Xi                       = xiAOD->Chi2Xi();
3199         lDcaXiDaughters                 = xiAOD->DcaXiDaughters();
3200         lDcaBachToPrimVertexXi          = xiAOD->DcaBachToPrimVertex();
3201         lXiCosineOfPointingAngle        = xiAOD->CosPointingAngleXi( lBestPrimaryVtxPos[0],
3202                                                                      lBestPrimaryVtxPos[1],
3203                                                                      lBestPrimaryVtxPos[2] );
3204
3205         lPosXi[0] = xiAOD->DecayVertexXiX();
3206         lPosXi[1] = xiAOD->DecayVertexXiY();
3207         lPosXi[2] = xiAOD->DecayVertexXiZ();
3208
3209                                         // This value shouldn't change, whatever the working hyp. is : Xi-, Xi+, Omega-, Omega+
3210         lDcaV0DaughtersXi               = xiAOD->DcaV0Daughters();
3211         // lV0Chi2Xi                    = xiAOD->GetChi2V0();
3212         lV0toXiCosineOfPointingAngle    = xiAOD->CosPointingAngle( lPosXi );
3213
3214         lV0CosineOfPointingAngle        = xiAOD->CosPointingAngle( lBestPrimaryVtxPos );
3215
3216         lPosV0Xi[0] = xiAOD->DecayVertexV0X();
3217         lPosV0Xi[1] = xiAOD->DecayVertexV0Y();
3218         lPosV0Xi[2] = xiAOD->DecayVertexV0Z();
3219
3220
3221         lDcaV0ToPrimVertexXi            = xiAOD->DcaV0ToPrimVertex();
3222
3223         lDcaPosToPrimVertexXi           = xiAOD->DcaPosToPrimVertex();
3224
3225         lDcaNegToPrimVertexXi           = xiAOD->DcaNegToPrimVertex();
3226
3227         lChargeXi = xiAOD->ChargeXi();
3228
3229         // - Step II.3 : PID info
3230
3231         //-------------
3232
3233
3234         // 3.1 - PID Information
3235
3236
3237         // Combined VO-positive-daughter PID
3238
3239         // Combined bachelor PID
3240 /*        AliPID bachPidXi;       bachPidXi.SetPriors(    lPriorsGuessXi    );
3241         AliPID bachPidOmega;    bachPidOmega.SetPriors( lPriorsGuessOmega );
3242
3243         if ( bachTrackXi->IsOn(AliESDtrack::kESDpid) ) {  // Combined PID exists
3244                 Double_t r[10] = {0.}; bachTrackXi->GetESDpid(r);
3245                 bachPidXi.SetProbabilities(r);
3246                 bachPidOmega.SetProbabilities(r);
3247                 // Check if the bachelor track is a pion
3248                     ppionBach = bachPidXi.GetProbability(AliPID::kPion);
3249                 if (ppionBach > bachPidXi.GetProbability(AliPID::kElectron) &&
3250                     ppionBach > bachPidXi.GetProbability(AliPID::kMuon)     &&
3251                     ppionBach > bachPidXi.GetProbability(AliPID::kKaon)     &&
3252                     ppionBach > bachPidXi.GetProbability(AliPID::kProton)   )     lIsBachelorPion = kTRUE;
3253                 // Check if the bachelor track is a kaon
3254                     pkaonBach = bachPidOmega.GetProbability(AliPID::kKaon);
3255                 if (pkaonBach > bachPidOmega.GetProbability(AliPID::kElectron) &&
3256                     pkaonBach > bachPidOmega.GetProbability(AliPID::kMuon)     &&
3257                     pkaonBach > bachPidOmega.GetProbability(AliPID::kPion)     &&
3258                     pkaonBach > bachPidOmega.GetProbability(AliPID::kProton)   )  lIsBachelorKaon = kTRUE;
3259         }// end if bachelor track with existing combined PID
3260 */
3261
3262         // 3.1.B - TPC PID : 4-sigma bands on Bethe-Bloch curve
3263
3264
3265                 // Bachelor
3266         if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( bachTrackXi,AliPID::kKaon)) < 4) lIsBachelorKaonForTPC = kTRUE;
3267         if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( bachTrackXi,AliPID::kPion)) < 4) lIsBachelorPionForTPC = kTRUE;
3268
3269         // Negative V0 daughter
3270         if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( nTrackXi,AliPID::kPion   )) < 4) lIsNegPionForTPC   = kTRUE;
3271         if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( nTrackXi,AliPID::kProton )) < 4) lIsNegProtonForTPC = kTRUE;
3272
3273         // Positive V0 daughter
3274         if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( pTrackXi,AliPID::kPion   )) < 4) lIsPosPionForTPC   = kTRUE;
3275         if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( pTrackXi,AliPID::kProton )) < 4) lIsPosProtonForTPC = kTRUE;
3276
3277 /*
3278         const AliExternalTrackParam *pInnerWallTrackXi    = pTrackXi    ->GetInnerParam(); // Do not use GetTPCInnerWall
3279         const AliExternalTrackParam *nInnerWallTrackXi    = nTrackXi    ->GetInnerParam();
3280         const AliExternalTrackParam *bachInnerWallTrackXi = bachTrackXi ->GetInnerParam();
3281         if(pInnerWallTrackXi && nInnerWallTrackXi && bachInnerWallTrackXi ){
3282
3283                 Double_t pMomInnerWall    = pInnerWallTrackXi   ->GetP();
3284                 Double_t nMomInnerWall    = nInnerWallTrackXi   ->GetP();
3285                 Double_t bachMomInnerWall = bachInnerWallTrackXi->GetP();
3286
3287                 // Bachelor
3288                 if (TMath::Abs(fESDpid->NumberOfSigmasTPC( bachTrackXi,AliPID::kPion)) < 3)                              lIsBachelorPionForTPC = kTRUE;
3289                 if (bachMomInnerWall < 0.350  && TMath::Abs(fESDpid->NumberOfSigmasTPC( bachTrackXi,AliPID::kKaon)) < 5) lIsBachelorKaonForTPC = kTRUE;
3290                 if (bachMomInnerWall > 0.350  && TMath::Abs(fESDpid->NumberOfSigmasTPC( bachTrackXi,AliPID::kKaon)) < 3) lIsBachelorKaonForTPC = kTRUE;
3291
3292                 // Negative V0 daughter
3293                 if (TMath::Abs(fESDpid->NumberOfSigmasTPC( nTrackXi,AliPID::kPion   )) < 3  )                           lIsNegPionForTPC   = kTRUE;
3294                 if (nMomInnerWall < 0.6  && TMath::Abs(fESDpid->NumberOfSigmasTPC( nTrackXi,AliPID::kProton ) ) < 5 )   lIsNegProtonForTPC = kTRUE;
3295                 if (nMomInnerWall > 0.6  && TMath::Abs(fESDpid->NumberOfSigmasTPC( nTrackXi,AliPID::kProton ) ) < 3 )   lIsNegProtonForTPC = kTRUE;
3296
3297                 // Positive V0 daughter
3298                 if (TMath::Abs(fESDpid->NumberOfSigmasTPC( pTrackXi,AliPID::kPion   )) < 3 )                            lIsPosPionForTPC   = kTRUE;
3299                 if (pMomInnerWall < 0.6  && TMath::Abs(fESDpid->NumberOfSigmasTPC( pTrackXi,AliPID::kProton )) < 5)     lIsPosProtonForTPC = kTRUE;
3300                 if (pMomInnerWall > 0.6  && TMath::Abs(fESDpid->NumberOfSigmasTPC( pTrackXi,AliPID::kProton )) < 3)     lIsPosProtonForTPC = kTRUE;
3301
3302         }
3303 */
3304
3305         // 3.2 - PID proba Vs Pt(Bach)
3306         Int_t      lblBachForPID  = (Int_t) TMath::Abs( bachTrackXi->GetLabel() );
3307         AliAODMCParticle* mcBachForPID   = (AliAODMCParticle*) arrayMC->At( lblBachForPID );
3308         lmcPtBach      = mcBachForPID->Pt();
3309
3310         // 3.3 - MC perfect PID
3311
3312         if( mcBachForPID->PdgCode() == -211) lIsBachelorMCPiMinus = kTRUE;
3313         if( mcBachForPID->PdgCode() ==  211) lIsBachelorMCPiPlus  = kTRUE;
3314         if( mcBachForPID->PdgCode() == -321) lIsBachelorMCKMinus  = kTRUE;
3315         if( mcBachForPID->PdgCode() ==  321) lIsBachelorMCKPlus   = kTRUE;
3316
3317
3318
3319         // - Step II.4 : MC association (care : lots of "continue;" below this line)
3320         //-------------
3321
3322         if(fDebug > 5)
3323                 cout    << "MC EventNumber : " << lMCevent->Header()->GetEvent()
3324                         << " / MC event Number in Run : " << lMCevent->Header()->GetEventNrInRun() << endl;
3325
3326         // - Step 4.1 : level of the V0 daughters
3327
3328         Int_t lblPosV0Dghter = (Int_t) TMath::Abs( pTrackXi->GetLabel() );
3329                 // Abs value = needed ! question of quality track association ... (negative label when at least one cluster in the track is from a different particle)
3330         Int_t lblNegV0Dghter = (Int_t) TMath::Abs( nTrackXi->GetLabel() );
3331         AliAODMCParticle* mcPosV0Dghter = (AliAODMCParticle*) arrayMC->At( lblPosV0Dghter );
3332         AliAODMCParticle* mcNegV0Dghter = (AliAODMCParticle*) arrayMC->At( lblNegV0Dghter );
3333
3334
3335         // - Step 4.2 : level of the Xi daughters
3336
3337         Int_t lblMotherPosV0Dghter = mcPosV0Dghter->GetMother();
3338         Int_t lblMotherNegV0Dghter = mcNegV0Dghter->GetMother();
3339
3340                 if( lblMotherPosV0Dghter != lblMotherNegV0Dghter) continue; // same mother
3341                 if( lblMotherPosV0Dghter < 0 ) continue; // this particle is primary, no mother
3342                 if( lblMotherNegV0Dghter < 0 ) continue;
3343
3344
3345                 // mothers = Lambda candidate ... a priori
3346
3347         AliAODMCParticle* mcMotherPosV0Dghter = (AliAODMCParticle*) arrayMC->At( lblMotherPosV0Dghter );
3348         AliAODMCParticle* mcMotherNegV0Dghter = (AliAODMCParticle*) arrayMC->At( lblMotherNegV0Dghter );  // MN: redundant?? already checked that labels are the same...-->same part from stack
3349
3350         Int_t      lblBach  = (Int_t) TMath::Abs( bachTrackXi->GetLabel() );
3351         AliAODMCParticle* mcBach   = (AliAODMCParticle*) arrayMC->At( lblBach );
3352
3353
3354         // - Step 4.3 : level of Xi candidate
3355
3356         Int_t lblGdMotherPosV0Dghter =   mcMotherPosV0Dghter->GetMother() ;
3357         Int_t lblGdMotherNegV0Dghter =   mcMotherNegV0Dghter->GetMother() ;
3358
3359                 if( lblGdMotherPosV0Dghter != lblGdMotherNegV0Dghter ) continue;
3360                 if( lblGdMotherPosV0Dghter < 0 ) continue; // primary lambda ...
3361                 if( lblGdMotherNegV0Dghter < 0 ) continue; // primary lambda ...
3362
3363                 // Gd mothers = Xi candidate ... a priori
3364
3365         AliAODMCParticle* mcGdMotherPosV0Dghter = (AliAODMCParticle*) arrayMC->At( lblGdMotherPosV0Dghter );
3366         AliAODMCParticle* mcGdMotherNegV0Dghter = (AliAODMCParticle*) arrayMC->At( lblGdMotherNegV0Dghter );
3367
3368         Int_t lblMotherBach = (Int_t) TMath::Abs( mcBach->GetMother() );
3369
3370         if( lblMotherBach != lblGdMotherPosV0Dghter ) continue; //same mother for bach and V0 daughters
3371
3372         AliAODMCParticle* mcMotherBach = (AliAODMCParticle*) arrayMC->At( lblMotherBach );
3373
3374                 // Check if cascade is primary
3375
3376         if (!(mcMotherBach->IsPhysicalPrimary())) continue;
3377
3378
3379         kIsNaturalPart = kTRUE;
3380
3381         if ((lblMotherBach>=endOfHijingEvent)&&(endOfHijingEvent!=-1)&&(mcMotherBach->GetMother()<0)) kIsNaturalPart = kFALSE;
3382
3383         // - Step 4.4 : Manage boolean for association
3384
3385         if( mcMotherBach                ->GetPdgCode() ==   3312 &&
3386             mcGdMotherPosV0Dghter       ->GetPdgCode() ==   3312 &&
3387             mcGdMotherNegV0Dghter       ->GetPdgCode() ==   3312)    {  lAssoXiMinus = kTRUE;
3388                                                                         lInvMassLambdaAsCascDghter = xiAOD->MassLambda();
3389                                                                         cascadeMass = 1.321;
3390                                                                       //  Printf("Xi- asso current index = %n ", lblGdMotherPosV0Dghter);
3391                                                                         nAssoXiMinus++; }
3392
3393         else if( mcMotherBach           ->GetPdgCode() ==  -3312 &&
3394             mcGdMotherPosV0Dghter       ->GetPdgCode() ==  -3312 &&
3395             mcGdMotherNegV0Dghter       ->GetPdgCode() ==  -3312)    {  lAssoXiPlus = kTRUE;
3396                                                                         lInvMassLambdaAsCascDghter = xiAOD->MassAntiLambda();
3397                                                                         cascadeMass = 1.321;
3398                                                                         nAssoXiPlus++; }
3399
3400         else if( mcMotherBach           ->GetPdgCode() ==   3334 &&
3401             mcGdMotherPosV0Dghter       ->GetPdgCode() ==   3334 &&
3402             mcGdMotherNegV0Dghter       ->GetPdgCode() ==   3334)    {  lAssoOmegaMinus = kTRUE;
3403                                                                         lInvMassLambdaAsCascDghter = xiAOD->MassLambda();
3404                                                                         cascadeMass = 1.672;
3405                                                                         nAssoOmegaMinus++; }
3406
3407         else if( mcMotherBach           ->GetPdgCode() ==  -3334 &&
3408             mcGdMotherPosV0Dghter       ->GetPdgCode() ==  -3334 &&
3409             mcGdMotherNegV0Dghter       ->GetPdgCode() ==  -3334)    {  lAssoOmegaPlus = kTRUE;
3410                                                                         lInvMassLambdaAsCascDghter = xiAOD->MassAntiLambda();
3411                                                                         cascadeMass = 1.672;
3412                                                                         nAssoOmegaPlus++; }
3413         //-------------
3414
3415         lmcPt             = mcMotherBach->Pt();
3416         lmcRapCasc        = 0.5*TMath::Log( (mcMotherBach->E() + mcMotherBach->Pz()) /
3417                                                      (mcMotherBach->E() - mcMotherBach->Pz() +1.e-13) );
3418         lmcEta            = mcMotherBach->Eta();
3419         Float_t decayCascX = mcBach->Xv();
3420         Float_t decayCascY = mcBach->Yv();
3421         lmcTransvRadius   = TMath::Sqrt(decayCascX*decayCascX+decayCascY*decayCascY); // decay point of Xi, = the production vertex of Bachelor ...
3422
3423         TVector3 lmcTVect3Mom( mcMotherBach->Px(), mcMotherBach->Py(), mcMotherBach->Pz() );
3424
3425         Double_t xiMomX = xiAOD->MomXiX();
3426         Double_t xiMomY = xiAOD->MomXiY();
3427         Double_t xiMomZ = xiAOD->MomXiZ();
3428
3429         lrecoPt           = TMath::Sqrt( xiMomX*xiMomX   + xiMomY*xiMomY );
3430         lrecoTransvRadius = TMath::Sqrt( xiAOD->DecayVertexXiX() * xiAOD->DecayVertexXiX() + xiAOD->DecayVertexXiY() * xiAOD->DecayVertexXiY() );
3431
3432         TVector3 lrecoTVect3Mom( xiMomX, xiMomY, xiMomZ );
3433         lDeltaPhiMcReco   = lmcTVect3Mom.DeltaPhi( lrecoTVect3Mom ) * 180.0/TMath::Pi();
3434
3435         lmcPtPosV0Dghter = mcPosV0Dghter->Pt() ;
3436         lmcPtNegV0Dghter = mcNegV0Dghter->Pt();
3437         lrecoP           = TMath::Sqrt( xiMomX*xiMomX   + xiMomY*xiMomY   + xiMomZ*xiMomZ );;
3438
3439         Double_t lV0momX = xiAOD->MomV0X();
3440         Double_t lV0momY = xiAOD->MomV0Y();
3441         Double_t lV0momZ = xiAOD->MomV0Z();
3442         lV0mom = TMath::Sqrt(TMath::Power(lV0momX,2)+TMath::Power(lV0momY,2)+TMath::Power(lV0momZ,2));
3443
3444       }
3445
3446       lXiRadius                       = TMath::Sqrt( lPosXi[0]*lPosXi[0]  +  lPosXi[1]*lPosXi[1] );
3447       lV0RadiusXi                     = TMath::Sqrt( lPosV0Xi[0]*lPosV0Xi[0]  +  lPosV0Xi[1]*lPosV0Xi[1] );
3448
3449       // Extra-selection for cascade candidates
3450                 // Towards optimisation of AA selection
3451       if (fkExtraSelections) {
3452
3453         // if(lChi2Xi > 2000) continue; // in AliCascadeVertexer
3454         // if(lV0Chi2Xi > 2000) continue; // in AliV0vertexer
3455
3456         if (lDcaXiDaughters > 0.3) continue; // in AliCascadeVertexer
3457         if (lXiCosineOfPointingAngle < 0.999 ) continue; // in AliCascadeVertexer
3458         if (lDcaV0ToPrimVertexXi < 0.05) continue; // in AliCascadeVertexer
3459         if (lDcaBachToPrimVertexXi < 0.03) continue; // in AliCascadeVertexer
3460 ////      if (TMath::Abs(lInvMassLambdaAsCascDghter-1.11568) > 0.006 ) continue;  // in AliCascadeVertexer
3461
3462         if (lDcaV0DaughtersXi > 1.) continue; // in AliV0vertexer
3463         if (lV0toXiCosineOfPointingAngle < 0.998) continue; // in AliV0vertexer
3464         if (lDcaPosToPrimVertexXi < 0.1) continue; // in AliV0vertexer
3465         if (lDcaNegToPrimVertexXi < 0.1) continue; // in AliV0vertexer
3466
3467
3468           if(lXiRadius < .9) continue; // in AliCascadeVertexer
3469 //        if(lXiRadius > 100) continue; // in AliCascadeVertexer
3470           if(lV0RadiusXi < 0.9) continue; // in AliV0vertexer
3471 //        if(lV0RadiusXi > 100) continue; // in AliV0vertexer
3472
3473       }
3474
3475
3476       // Combined PID TH1s
3477       if( lChargeXi < 0 && lIsBachelorPion )    fHistMassWithCombPIDXiMinus     ->Fill( lInvMassXiMinus    );
3478       if( lChargeXi > 0 && lIsBachelorPion )    fHistMassWithCombPIDXiPlus      ->Fill( lInvMassXiPlus     );
3479       if( lChargeXi < 0 && lIsBachelorKaon )    fHistMassWithCombPIDOmegaMinus  ->Fill( lInvMassOmegaMinus );
3480       if( lChargeXi > 0 && lIsBachelorKaon )    fHistMassWithCombPIDOmegaPlus   ->Fill( lInvMassOmegaPlus  );
3481
3482
3483       if( lChargeXi < 0 )     fHistMassXiMinus        ->Fill( lInvMassXiMinus );
3484       if( lChargeXi > 0 )     fHistMassXiPlus         ->Fill( lInvMassXiPlus );
3485       if( lChargeXi < 0 )     fHistMassOmegaMinus     ->Fill( lInvMassOmegaMinus );
3486       if( lChargeXi > 0 )     fHistMassOmegaPlus      ->Fill( lInvMassOmegaPlus );
3487
3488       if(lIsBachelorPion)   f2dHistPIDprobaPionVsMCPtBach->Fill( lmcPtBach, ppionBach );
3489       if(lIsBachelorKaon)   f2dHistPIDprobaKaonVsMCPtBach->Fill( lmcPtBach, pkaonBach );
3490
3491       if( lChargeXi < 0 && lIsBachelorMCPiMinus )    fHistMassWithMcPIDXiMinus     ->Fill( lInvMassXiMinus );
3492       if( lChargeXi > 0 && lIsBachelorMCPiPlus  )    fHistMassWithMcPIDXiPlus      ->Fill( lInvMassXiPlus );
3493       if( lChargeXi < 0 && lIsBachelorMCKMinus  )    fHistMassWithMcPIDOmegaMinus  ->Fill( lInvMassOmegaMinus );
3494       if( lChargeXi > 0 && lIsBachelorMCKPlus   )    fHistMassWithMcPIDOmegaPlus   ->Fill( lInvMassOmegaPlus );
3495
3496
3497
3498       if(!lAssoXiMinus && !lAssoXiPlus && !lAssoOmegaMinus && !lAssoOmegaPlus) continue; // no association, skip the rest of the code
3499
3500         // Calculate proper time for cascade (reconstructed)
3501
3502         Double_t lctau =  TMath::Sqrt(TMath::Power((lPosXi[0]-lBestPrimaryVtxPos[0]),2)+TMath::Power((lPosXi[1]-lBestPrimaryVtxPos[1]),2)+TMath::Power((lPosXi[2]-lBestPrimaryVtxPos[2]),2));
3503         if (lrecoP!=0)         lctau = lctau*cascadeMass/lrecoP;
3504         else lctau = -1.;
3505
3506
3507         // Calculate proper time for Lambda (reconstructed)
3508         Float_t lambdaMass = 1.115683; // PDG mass
3509         Float_t distV0Xi =  TMath::Sqrt(TMath::Power((lPosV0Xi[0]-lPosXi[0]),2)+TMath::Power((lPosV0Xi[1]-lPosXi[1]),2)+TMath::Power((lPosV0Xi[2]-lPosXi[2]),2));
3510         Float_t lctauV0 = -1.;
3511         if (lV0mom!=0) lctauV0 = distV0Xi*lambdaMass/lV0mom;
3512
3513         Float_t distTV0Xi =  TMath::Sqrt(TMath::Power((lPosV0Xi[0]-lPosXi[0]),2)+TMath::Power((lPosV0Xi[1]-lPosXi[1]),2));
3514
3515
3516         // - Histos for the cascade candidates associated with MC
3517
3518         if( lChargeXi < 0 && lAssoXiMinus){
3519                 fHistAsMCMassXiMinus          ->Fill( lInvMassXiMinus  );
3520                 if(lIsBachelorPion)     f2dHistAsMCandCombPIDGenPtVsGenYXiMinus->Fill( lmcPt, lmcRapCasc );
3521                 f2dHistAsMCGenPtVsGenYXiMinus ->Fill( lmcPt, lmcRapCasc);
3522                 fHistAsMCGenEtaXiMinus        ->Fill( lmcEta           );
3523                 f2dHistAsMCResPtXiMinus       ->Fill( lmcPt,           (lrecoPt - lmcPt)/ lmcPt );
3524                 f2dHistAsMCResRXiMinus        ->Fill( lmcTransvRadius, (lrecoTransvRadius - lmcTransvRadius)/ lmcTransvRadius    );
3525                 f2dHistAsMCResPhiXiMinus      ->Fill( lmcPt, lDeltaPhiMcReco );
3526                 f2dHistAsMCptProtonMCptXiMinus->Fill(lmcPt,lmcPtPosV0Dghter);
3527                 fHistV0CosineOfPointingAnglevsPtXi->Fill(lmcPt,lV0CosineOfPointingAngle);
3528         }
3529
3530         else if( lChargeXi > 0 && lAssoXiPlus){
3531                 fHistAsMCMassXiPlus           ->Fill( lInvMassXiPlus   );
3532                 if(lIsBachelorPion)     f2dHistAsMCandCombPIDGenPtVsGenYXiPlus->Fill( lmcPt, lmcRapCasc );
3533                 f2dHistAsMCGenPtVsGenYXiPlus  ->Fill( lmcPt, lmcRapCasc);
3534                 fHistAsMCGenEtaXiPlus         ->Fill( lmcEta           );
3535                 f2dHistAsMCResPtXiPlus        ->Fill( lmcPt,           (lrecoPt - lmcPt)/ lmcPt );
3536                 f2dHistAsMCResRXiPlus         ->Fill( lmcTransvRadius, (lrecoTransvRadius - lmcTransvRadius)/ lmcTransvRadius    );
3537                 f2dHistAsMCResPhiXiPlus       ->Fill( lmcPt, lDeltaPhiMcReco );
3538                 f2dHistAsMCptAntiprotonMCptXiPlus->Fill(lmcPt,lmcPtNegV0Dghter);
3539                 fHistV0CosineOfPointingAnglevsPtXi->Fill(lmcPt,lV0CosineOfPointingAngle);
3540         }
3541
3542         else if( lChargeXi < 0 && lAssoOmegaMinus){
3543                 fHistAsMCMassOmegaMinus          ->Fill( lInvMassOmegaMinus );
3544                 if(lIsBachelorKaon)     f2dHistAsMCandCombPIDGenPtVsGenYOmegaMinus->Fill( lmcPt, lmcRapCasc );
3545                 f2dHistAsMCGenPtVsGenYOmegaMinus ->Fill( lmcPt, lmcRapCasc  );
3546                 fHistAsMCGenEtaOmegaMinus        ->Fill( lmcEta             );
3547                 f2dHistAsMCResPtOmegaMinus       ->Fill( lmcPt,           (lrecoPt - lmcPt)/ lmcPt );
3548                 f2dHistAsMCResROmegaMinus        ->Fill( lmcTransvRadius, (lrecoTransvRadius - lmcTransvRadius)/ lmcTransvRadius    );
3549                 f2dHistAsMCResPhiOmegaMinus      ->Fill( lmcPt, lDeltaPhiMcReco );
3550                 f2dHistAsMCptProtonMCptOmegaMinus->Fill(lmcPt,lmcPtPosV0Dghter);
3551                 fHistV0CosineOfPointingAnglevsPtOmega->Fill(lmcPt,lV0CosineOfPointingAngle);
3552         }
3553
3554         else if( lChargeXi > 0 && lAssoOmegaPlus){
3555                 fHistAsMCMassOmegaPlus           ->Fill( lInvMassOmegaPlus );
3556                 if(lIsBachelorKaon)     f2dHistAsMCandCombPIDGenPtVsGenYOmegaPlus->Fill( lmcPt, lmcRapCasc );
3557                 f2dHistAsMCGenPtVsGenYOmegaPlus  ->Fill( lmcPt, lmcRapCasc   );
3558                 fHistAsMCGenEtaOmegaPlus         ->Fill( lmcEta            );
3559                 f2dHistAsMCResPtOmegaPlus        ->Fill( lmcPt,           (lrecoPt - lmcPt)/ lmcPt );
3560                 f2dHistAsMCResROmegaPlus         ->Fill( lmcTransvRadius, (lrecoTransvRadius - lmcTransvRadius)/ lmcTransvRadius    );
3561                 f2dHistAsMCResPhiOmegaPlus       ->Fill( lmcPt, lDeltaPhiMcReco );
3562                 f2dHistAsMCptAntiprotonMCptOmegaPlus->Fill(lmcPt,lmcPtNegV0Dghter);
3563                 fHistV0CosineOfPointingAnglevsPtOmega->Fill(lmcPt,lV0CosineOfPointingAngle);
3564         }
3565
3566         fHistV0toXiCosineOfPointingAngle->Fill(lV0toXiCosineOfPointingAngle);
3567
3568         // - Step 6 : Containers = Cascade cuts + PID
3569         //-------------
3570         // 6.3 - Filling the AliCFContainer (optimisation of topological selections + systematics)
3571         Double_t lContainerCutVars[21] = {0.0};
3572
3573         lContainerCutVars[0]  = lDcaXiDaughters;
3574         lContainerCutVars[1]  = lDcaBachToPrimVertexXi;
3575         lContainerCutVars[2]  = lXiCosineOfPointingAngle;
3576         lContainerCutVars[3]  = lXiRadius;
3577         lContainerCutVars[4]  = lInvMassLambdaAsCascDghter;
3578         lContainerCutVars[5]  = lDcaV0DaughtersXi;
3579         lContainerCutVars[6]  = lV0CosineOfPointingAngle;
3580         lContainerCutVars[7]  = lV0RadiusXi;
3581         lContainerCutVars[8]  = lDcaV0ToPrimVertexXi;
3582         lContainerCutVars[9]  = lDcaPosToPrimVertexXi;
3583         lContainerCutVars[10] = lDcaNegToPrimVertexXi;
3584
3585         lContainerCutVars[13] = lmcPt;
3586
3587         if (kIsNaturalPart) lContainerCutVars[16] = 0.;
3588         else lContainerCutVars[16] = 1.;
3589         lContainerCutVars[17] = lcentrality;
3590         lContainerCutVars[18] = lPrimaryTrackMultiplicity;
3591         lContainerCutVars[19] = lctau;
3592         lContainerCutVars[20] = lctauV0;
3593         lContainerCutVars[21] = distTV0Xi;
3594
3595         // All cases should be covered below
3596         if( lChargeXi < 0 && lAssoXiMinus    ) {
3597                 lContainerCutVars[11] = lInvMassXiMinus;
3598                 lContainerCutVars[12] = lInvMassOmegaMinus;//1.63;
3599                 lContainerCutVars[14] = lmcRapCasc;
3600                 lContainerCutVars[15] = -1.;
3601                 if ( lIsBachelorPionForTPC   && lIsPosProtonForTPC    && lIsNegPionForTPC )
3602                   fCFContAsCascadeCuts->Fill(lContainerCutVars,0); // for Xi-
3603         }
3604         if( lChargeXi > 0 && lAssoXiPlus    ) {
3605                 lContainerCutVars[11] = lInvMassXiPlus;
3606                 lContainerCutVars[12] = lInvMassOmegaPlus;//1.26;
3607                 lContainerCutVars[14] = lmcRapCasc;
3608                 lContainerCutVars[15] = -1.;
3609                 if ( lIsBachelorPionForTPC   && lIsNegProtonForTPC    && lIsPosPionForTPC )
3610                   fCFContAsCascadeCuts->Fill(lContainerCutVars,1); // for Xi+
3611         }
3612         if( lChargeXi < 0 && lAssoOmegaMinus ) {
3613                 lContainerCutVars[11] = lInvMassXiMinus;//1.63;
3614                 lContainerCutVars[12] = lInvMassOmegaMinus;
3615                 lContainerCutVars[14] = -1.;
3616                 lContainerCutVars[15] = lmcRapCasc;
3617                 if ( lIsBachelorKaonForTPC   && lIsPosProtonForTPC    && lIsNegPionForTPC )
3618                   fCFContAsCascadeCuts->Fill(lContainerCutVars,2); // for Omega-
3619         }
3620         if( lChargeXi > 0 && lAssoOmegaPlus  ) {
3621                 lContainerCutVars[11] = lInvMassXiPlus;//1.26;
3622                 lContainerCutVars[12] = lInvMassOmegaPlus;
3623                 lContainerCutVars[14] = -1.;
3624                 lContainerCutVars[15] = lmcRapCasc;
3625                 if ( lIsBachelorKaonForTPC   && lIsNegProtonForTPC    && lIsPosPionForTPC )
3626                   fCFContAsCascadeCuts->Fill(lContainerCutVars,3); // for Omega+
3627         }
3628
3629
3630         // 6.4 - Filling the AliCFContainers related to PID
3631
3632         Double_t lContainerPIDVars[4] = {0.0};
3633
3634
3635         // Xi Minus
3636         if( lChargeXi < 0 && lAssoXiMinus ) {
3637                 lContainerPIDVars[0] = lmcPt              ;
3638                 lContainerPIDVars[1] = lInvMassXiMinus    ;
3639                 lContainerPIDVars[2] = lmcRapCasc         ;
3640                 lContainerPIDVars[3] = lcentrality;
3641
3642                 // No PID
3643                         fCFContCascadePIDAsXiMinus->Fill(lContainerPIDVars, 0); // No PID
3644                 // TPC PID
3645                 if( lIsBachelorPionForTPC  )
3646                         fCFContCascadePIDAsXiMinus->Fill(lContainerPIDVars, 1); // TPC PID / 4-#sigma cut on Bachelor track
3647
3648                 if( lIsBachelorPionForTPC &&
3649                     lIsPosProtonForTPC     )
3650                         fCFContCascadePIDAsXiMinus->Fill(lContainerPIDVars, 2); // TPC PID / 4-#sigma cut on Bachelor+Baryon tracks
3651
3652                 if( lIsBachelorPionForTPC &&
3653                     lIsPosProtonForTPC    &&
3654                     lIsNegPionForTPC       )
3655                         fCFContCascadePIDAsXiMinus->Fill(lContainerPIDVars, 3); // TPC PID / 4-#sigma cut on Bachelor+Baryon+Meson tracks
3656
3657                 // Combined PID
3658                 if( lIsBachelorPion        )
3659                         fCFContCascadePIDAsXiMinus->Fill(lContainerPIDVars, 4); // Comb. PID / Bachelor
3660
3661                 if( lIsBachelorPion       &&
3662                     lIsPosInXiProton    )
3663                         fCFContCascadePIDAsXiMinus->Fill(lContainerPIDVars, 5); // Comb. PID / Bachelor+Baryon
3664
3665                 if( lIsBachelorPion     &&
3666                     lIsPosInXiProton &&
3667                     lIsNegInXiPion    )
3668                         fCFContCascadePIDAsXiMinus->Fill(lContainerPIDVars, 6); // Comb. PID / Bachelor+Baryon+Meson
3669         }
3670
3671         lContainerPIDVars[0] = 0.; lContainerPIDVars[1] = 0.; lContainerPIDVars[2] = 0.; lContainerPIDVars[3] = 0.;
3672
3673         // Xi Plus
3674         if( lChargeXi > 0 && lAssoXiPlus ) {
3675                 lContainerPIDVars[0] = lmcPt         ;
3676                 lContainerPIDVars[1] = lInvMassXiPlus;
3677                 lContainerPIDVars[2] = lmcRapCasc    ;
3678                 lContainerPIDVars[3] = lcentrality   ;
3679
3680                 // No PID
3681                         fCFContCascadePIDAsXiPlus->Fill(lContainerPIDVars, 0); // No PID
3682                 // TPC PID
3683                 if( lIsBachelorPionForTPC  )
3684                         fCFContCascadePIDAsXiPlus->Fill(lContainerPIDVars, 1); // TPC PID / 4-#sigma cut on Bachelor track
3685
3686                 if( lIsBachelorPionForTPC &&
3687                     lIsNegProtonForTPC     )
3688                         fCFContCascadePIDAsXiPlus->Fill(lContainerPIDVars, 2); // TPC PID / 4-#sigma cut on Bachelor+Baryon tracks
3689
3690                 if( lIsBachelorPionForTPC &&
3691                     lIsNegProtonForTPC    &&
3692                     lIsPosPionForTPC       )
3693                         fCFContCascadePIDAsXiPlus->Fill(lContainerPIDVars, 3); // TPC PID / 4-#sigma cut on Bachelor+Baryon+Meson tracks
3694
3695                 // Combined PID
3696                 if( lIsBachelorPion        )
3697                         fCFContCascadePIDAsXiPlus->Fill(lContainerPIDVars, 4); // Comb. PID / Bachelor
3698
3699                 if( lIsBachelorPion       &&
3700                     lIsNegInXiProton    )
3701                         fCFContCascadePIDAsXiPlus->Fill(lContainerPIDVars, 5); // Comb. PID / Bachelor+Baryon
3702
3703                 if( lIsBachelorPion     &&
3704                     lIsNegInXiProton &&
3705                     lIsPosInXiPion    )
3706                         fCFContCascadePIDAsXiPlus->Fill(lContainerPIDVars, 6); // Comb. PID / Bachelor+Baryon+Meson
3707         }
3708
3709         lContainerPIDVars[0] = 0.; lContainerPIDVars[1] = 0.; lContainerPIDVars[2] = 0.; lContainerPIDVars[3] = 0.;
3710
3711         // Omega Minus
3712         if( lChargeXi < 0 && lAssoOmegaMinus ) {
3713                 lContainerPIDVars[0] = lmcPt              ;
3714                 lContainerPIDVars[1] = lInvMassOmegaMinus ;
3715                 lContainerPIDVars[2] = lmcRapCasc         ;
3716                 lContainerPIDVars[3] = lcentrality;
3717
3718                 // No PID
3719                         fCFContCascadePIDAsOmegaMinus->Fill(lContainerPIDVars, 0); // No PID
3720                 // TPC PID
3721                 if( lIsBachelorKaonForTPC  )
3722                         fCFContCascadePIDAsOmegaMinus->Fill(lContainerPIDVars, 1); // TPC PID / 4-#sigma cut on Bachelor track
3723
3724                 if( lIsBachelorKaonForTPC &&
3725                     lIsPosProtonForTPC     )
3726                         fCFContCascadePIDAsOmegaMinus->Fill(lContainerPIDVars, 2); // TPC PID / 4-#sigma cut on Bachelor+Baryon tracks
3727
3728                 if( lIsBachelorKaonForTPC &&
3729                     lIsPosProtonForTPC    &&
3730                     lIsNegPionForTPC       )
3731                         fCFContCascadePIDAsOmegaMinus->Fill(lContainerPIDVars, 3); // TPC PID / 4-#sigma cut on Bachelor+Baryon+Meson tracks
3732
3733                 // Combined PID
3734                 if( lIsBachelorKaon        )
3735                         fCFContCascadePIDAsOmegaMinus->Fill(lContainerPIDVars, 4); // Comb. PID / Bachelor
3736
3737                 if( lIsBachelorKaon       &&
3738                     lIsPosInOmegaProton    )
3739                         fCFContCascadePIDAsOmegaMinus->Fill(lContainerPIDVars, 5); // Comb. PID / Bachelor+Baryon
3740
3741                 if(lIsBachelorKaon     &&
3742                    lIsPosInOmegaProton &&
3743                    lIsNegInOmegaPion    )
3744                         fCFContCascadePIDAsOmegaMinus->Fill(lContainerPIDVars, 6); // Comb. PID / Bachelor+Baryon+Meson
3745         }
3746
3747         lContainerPIDVars[0] = 0.; lContainerPIDVars[1] = 0.; lContainerPIDVars[2] = 0.; lContainerPIDVars[3] = 0.;
3748
3749         // Omega Plus
3750         if( lChargeXi > 0 && lAssoOmegaPlus) {
3751                 lContainerPIDVars[0] = lmcPt              ;
3752                 lContainerPIDVars[1] = lInvMassOmegaPlus  ;
3753                 lContainerPIDVars[2] = lmcRapCasc         ;
3754                 lContainerPIDVars[3] = lcentrality;
3755
3756                 // No PID
3757                         fCFContCascadePIDAsOmegaPlus->Fill(lContainerPIDVars, 0); // No PID
3758                 // TPC PID
3759                 if( lIsBachelorKaonForTPC  )
3760                         fCFContCascadePIDAsOmegaPlus->Fill(lContainerPIDVars, 1); // TPC PID / 4-#sigma cut on Bachelor track
3761
3762                 if( lIsBachelorKaonForTPC &&
3763                     lIsNegProtonForTPC     )
3764                         fCFContCascadePIDAsOmegaPlus->Fill(lContainerPIDVars, 2); // TPC PID / 4-#sigma cut on Bachelor+Baryon tracks
3765
3766                 if( lIsBachelorKaonForTPC &&
3767                     lIsNegProtonForTPC    &&
3768                     lIsPosPionForTPC       )
3769                         fCFContCascadePIDAsOmegaPlus->Fill(lContainerPIDVars, 3); // TPC PID / 4-#sigma cut on Bachelor+Baryon+Meson tracks
3770
3771                 // Combined PID
3772                 if( lIsBachelorKaon        )
3773                         fCFContCascadePIDAsOmegaPlus->Fill(lContainerPIDVars, 4); // Comb. PID / Bachelor
3774
3775                 if( lIsBachelorKaon       &&
3776                     lIsNegInOmegaProton    )
3777                         fCFContCascadePIDAsOmegaPlus->Fill(lContainerPIDVars, 5); // Comb. PID / Bachelor+Baryon
3778
3779                 if( lIsBachelorKaon     &&
3780                     lIsNegInOmegaProton &&
3781                     lIsPosInOmegaPion    )
3782                         fCFContCascadePIDAsOmegaPlus->Fill(lContainerPIDVars, 6); // Comb. PID / Bachelor+Baryon+Meson
3783
3784       }
3785
3786 }// End of loop over reconstructed cascades
3787
3788 fHistnAssoXiMinus->Fill(nAssoXiMinus);
3789 fHistnAssoXiPlus->Fill(nAssoXiPlus);
3790 fHistnAssoOmegaMinus->Fill(nAssoOmegaMinus);
3791 fHistnAssoOmegaPlus->Fill(nAssoOmegaPlus);
3792 //Printf("N asso Xi- = %n ", nAssoXiMinus);
3793
3794   // Post output data.
3795  PostData(1, fListHistCascade);
3796 }
3797
3798
3799 //________________________________________________________________________
3800 void AliAnalysisTaskCheckPerformanceCascadePbPb::Terminate(Option_t *) {
3801   // Draw result to the screen
3802   // Called once at the end of the query
3803
3804   TList *cRetrievedList = 0x0;
3805   cRetrievedList = (TList*)GetOutputData(1);
3806   if(!cRetrievedList) {
3807         Printf("ERROR - AliAnalysisTaskCheckPerformanceCascadePbPb : ouput data container list not available\n");
3808         return;
3809   }
3810
3811   fHistMCTrackMultiplicity = dynamic_cast<TH1F*> (  cRetrievedList->FindObject("fHistMCTrackMultiplicity")  );
3812   if (!fHistMCTrackMultiplicity) {
3813     Printf("ERROR - AliAnalysisTaskCheckPerformanceCascadePbPb : fHistMCTrackMultiplicity not available");
3814     return;
3815   }
3816
3817
3818   TCanvas *canCheckPerformanceCascade = new TCanvas("AliAnalysisTaskCheckPerformanceCascadePbPb","Multiplicity",10,10,510,510);
3819   canCheckPerformanceCascade->cd(1)->SetLogy();
3820
3821   fHistMCTrackMultiplicity->SetMarkerStyle(22);
3822   fHistMCTrackMultiplicity->DrawCopy("E");
3823
3824 }