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