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 = dynamic_cast<AliAODTrack*>(lAODevent->GetTrack(itrack));
1970             if(!track) AliFatal("Not a standard AOD");
1971             if (track->TestFilterBit(AliAODTrack::kTrkGlobalNoDCA)) lPrimaryTrackMultiplicity++;
1972           }
1973
1974 //          lSPDTrackletsMultiplicity = lAODevent->GetTracklets()->GetNumberOfTracklets();
1975
1976           nNumberOfMCPrimaries = arrayMC->GetEntries();
1977
1978
1979           centrality = lAODevent->GetCentrality();
1980           aodV0 = lAODevent->GetVZEROData();
1981           multV0A=aodV0->GetMTotV0A();
1982           multV0C=aodV0->GetMTotV0C();
1983
1984         }
1985
1986         if (fkQualityCutZprimVtxPos) { // Should be already in the centrality selection
1987           if (TMath::Abs(lBestPrimaryVtxPos[2]) > fVtxRange ) {
1988             AliWarning("Pb / | Z position of Best Prim Vtx | > 10.0 cm ... return !");
1989             PostData(1, fListHistCascade);
1990             PostData(2, fCFContCascadePIDAsXiMinus);
1991             PostData(3, fCFContCascadePIDAsXiPlus);
1992             PostData(4, fCFContCascadePIDAsOmegaMinus);
1993             PostData(5, fCFContCascadePIDAsOmegaPlus);
1994             PostData(6, fCFContAsCascadeCuts);
1995
1996             return;
1997           }
1998         }
1999
2000
2001       //  Printf("Centrality percentile V0M for this event %f)\n",  centrality->GetCentralityPercentile("V0M"));
2002
2003
2004         Float_t lcentrality = 0.;
2005
2006         if (fkUseCleaning) lcentrality = centrality->GetCentralityPercentile(fCentrEstimator.Data());
2007         else {
2008           lcentrality = centrality->GetCentralityPercentileUnchecked(fCentrEstimator.Data());
2009           if (centrality->GetQuality()>1) {
2010             PostData(1, fListHistCascade);
2011             PostData(2, fCFContCascadePIDAsXiMinus);
2012             PostData(3, fCFContCascadePIDAsXiPlus);
2013             PostData(4, fCFContCascadePIDAsOmegaMinus);
2014             PostData(5, fCFContCascadePIDAsOmegaPlus);
2015             PostData(6, fCFContAsCascadeCuts);
2016
2017             return;
2018           }
2019         }
2020         if (lcentrality<fCentrLowLim||lcentrality>=fCentrUpLim) {
2021           PostData(1, fListHistCascade);
2022           return;
2023         }
2024
2025         if (nNumberOfMCPrimaries < 1) {
2026           PostData(1, fListHistCascade);
2027           PostData(2, fCFContCascadePIDAsXiMinus);
2028           PostData(3, fCFContCascadePIDAsXiPlus);
2029           PostData(4, fCFContCascadePIDAsOmegaMinus);
2030           PostData(5, fCFContCascadePIDAsOmegaPlus);
2031           PostData(6, fCFContAsCascadeCuts);
2032
2033           return; // should be useless because we require vertex and centrality selection
2034         }
2035
2036         fV0Ampl->Fill(multV0A+multV0C);
2037
2038         // Best vertex distribution
2039         fHistBestVtxX->Fill(lBestPrimaryVtxPos[0]);
2040         fHistBestVtxY->Fill(lBestPrimaryVtxPos[1]);
2041         fHistBestVtxZ->Fill(lBestPrimaryVtxPos[2]);
2042
2043         fHistEvtsInCentralityBinsvsNtracks->Fill(lcentrality,lPrimaryTrackMultiplicity);
2044
2045         fHistMCTrackMultiplicity->Fill( nNumberOfMCPrimaries );  // MN: neutral particles included and also not physical ones
2046
2047         //cout << "Name of the accessed file :" <<  CurrentFileName() << endl;
2048
2049         Int_t   nMCPrimariesInAcceptance   =  0;
2050 //_____________________________________________________________________________
2051 // Part 1 - Loop over the MC primaries
2052
2053     for (Int_t iCurrentLabelStack = 0; iCurrentLabelStack < nNumberOfMCPrimaries; iCurrentLabelStack++) {
2054
2055       if (fAnalysisType == "ESD") {
2056
2057         TParticle* lCurrentParticle = 0x0;
2058                    lCurrentParticle = lMCstack->Particle( iCurrentLabelStack );
2059         if (!lCurrentParticle) {
2060           Printf("MC Primary loop %d - MC TParticle pointer to current stack particle = 0x0 ! Skip ...\n", iCurrentLabelStack );
2061           continue;
2062         }
2063
2064         if (!lMCstack->IsPhysicalPrimary(iCurrentLabelStack)) continue;
2065
2066
2067         if ( TMath::Abs( lCurrentParticle->Eta() ) > 0.8 ) continue;
2068
2069       } else if (fAnalysisType == "AOD") {
2070
2071         AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iCurrentLabelStack);
2072         if(TMath::Abs(partMC->Eta())>0.8) continue;
2073
2074       }
2075
2076       nMCPrimariesInAcceptance++;
2077
2078    }
2079
2080    f2dHistRecoMultVsMCMult->Fill( lPrimaryTrackMultiplicity, nMCPrimariesInAcceptance );  // MN: neutral are included
2081
2082
2083    // For proton
2084
2085 /*   for (Int_t iCurrentLabelStack = 0; iCurrentLabelStack < nNumberOfMCPrimaries; iCurrentLabelStack++) {
2086
2087         TParticle* lCurrentParticle = lMCstack->Particle( iCurrentLabelStack );
2088         if(!lCurrentParticle){
2089                 Printf("Proton loop %d - MC TParticle pointer to current stack particle = 0x0 ! Skip ...\n", iCurrentLabelStack );
2090                 continue;
2091
2092         }
2093
2094         if( lCurrentParticle->GetPdgCode() == 2212 )
2095                 fHistEtaGenProton->Fill( lCurrentParticle->Eta() );
2096
2097         if( lCurrentParticle->GetPdgCode() == -2212 )
2098                 fHistEtaGenAntiProton->Fill( lCurrentParticle->Eta() );
2099       }// end loop over primary proton
2100 */
2101
2102
2103
2104 //_____________________________________________________________________________
2105 // Part 2 - Loop over the different types of GENERATED cascades (Xi-+, Omega-+)
2106
2107         // - Initialisation of useful local variables
2108
2109         Int_t lPdgCodeCasc            = 0;
2110         Int_t lPdgCodeBach            = 0;
2111         Int_t lPdgCodeLambda          = 0;
2112         Int_t lPdgCodeDghtMesV0       = 0;
2113         Int_t lPdgCodeDghtBarV0       = 0;
2114
2115
2116         TH1F *lHistEtaGenCasc         = 0;
2117         TH3D *l3dHistGenPtVsGenYvsCentNat = 0;
2118         TH3D *l3dHistGenPtVsGenYvsNtracksNat = 0;
2119         TH3D *l3dHistGenPtVsGenYvsCentInj = 0;
2120         TH3D *l3dHistGenPtVsGenYvsNtracksInj = 0;
2121         TH3D *l3dHistGenPtVsGenctauvsCentNat = 0;
2122         TH3D *l3dHistGenPtVsGenctauvsCentInj = 0;
2123
2124         TH1F *lHistThetaGenCascNat    = 0;
2125         TH1F *lHistThetaGenCascInj    = 0;
2126         TH2D *l2dHistGenPtVsGenYFdbl  = 0;
2127         TH1F *lHistThetaLambda        = 0;
2128         TH1F *lHistThetaBach          = 0;
2129         TH1F *lHistThetaBarDghter     = 0;
2130         TH1F *lHistThetaMesDghter     = 0;
2131         TH1F *lHistPtBach             = 0;
2132         TH1F *lHistPtBarDghter        = 0;
2133         TH1F *lHistPtMesDghter        = 0;
2134         Int_t ncascperev = 0;
2135         Int_t ncascperevtot = 0;
2136
2137         Int_t endOfHijingEvent = -1;
2138
2139 for (Int_t iCascType = 1; iCascType < 5; iCascType++) {
2140   ncascperev = 0;
2141   ncascperevtot = 0;
2142
2143   endOfHijingEvent = -1;
2144
2145   switch (iCascType) {
2146     case 1: // Xi-
2147          lPdgCodeCasc       =   3312;  //Xi-
2148          lPdgCodeBach       =   -211;  //Pi-
2149          lPdgCodeLambda     =   3122;  //Lambda0
2150          lPdgCodeDghtMesV0  =   -211;  //Pi-
2151          lPdgCodeDghtBarV0  =   2212;  //Proton
2152
2153                 // any Xi-
2154          lHistEtaGenCasc        = fHistEtaGenCascXiMinus;
2155          l3dHistGenPtVsGenYvsCentNat = f3dHistGenPtVsGenYvsCentXiMinusNat;
2156          l3dHistGenPtVsGenYvsNtracksNat = f3dHistGenPtVsGenYvsNtracksXiMinusNat;
2157          l3dHistGenPtVsGenYvsCentInj = f3dHistGenPtVsGenYvsCentXiMinusInj;
2158          l3dHistGenPtVsGenYvsNtracksInj = f3dHistGenPtVsGenYvsNtracksXiMinusInj;
2159          l3dHistGenPtVsGenctauvsCentNat = f3dHistGenPtVsGenctauvsCentXiMinusNat;
2160          l3dHistGenPtVsGenctauvsCentInj = f3dHistGenPtVsGenctauvsCentXiMinusInj;
2161
2162                 // cascades generated within acceptance (cut in pt + theta)
2163          lHistThetaGenCascNat      = fHistThetaGenCascXiMinusNat;
2164          lHistThetaGenCascInj      = fHistThetaGenCascXiMinusInj;
2165          l2dHistGenPtVsGenYFdbl = f2dHistGenPtVsGenYFdblXiMinus;
2166          lHistThetaLambda       = fHistThetaLambdaXiMinus;
2167          lHistThetaBach         = fHistThetaBachXiMinus;
2168          lHistThetaBarDghter    = fHistThetaBarDghterXiMinus;
2169          lHistThetaMesDghter    = fHistThetaMesDghterXiMinus;
2170          lHistPtBach            = fHistPtBachXiMinus;
2171          lHistPtBarDghter       = fHistPtBarDghterXiMinus;
2172          lHistPtMesDghter       = fHistPtMesDghterXiMinus;
2173         break;
2174
2175     case 2: // Xi+
2176          lPdgCodeCasc        =  -3312;  //Xi+
2177          lPdgCodeBach        =    211;  //Pi+
2178          lPdgCodeLambda      =  -3122;  //AntiLambda0
2179          lPdgCodeDghtMesV0   =    211;  //Pi+
2180          lPdgCodeDghtBarV0   =  -2212;  //AntiProton
2181
2182                 // any Xi+
2183          lHistEtaGenCasc        = fHistEtaGenCascXiPlus;
2184          l3dHistGenPtVsGenYvsCentNat = f3dHistGenPtVsGenYvsCentXiPlusNat;
2185          l3dHistGenPtVsGenYvsNtracksNat = f3dHistGenPtVsGenYvsNtracksXiPlusNat;
2186          l3dHistGenPtVsGenYvsCentInj = f3dHistGenPtVsGenYvsCentXiPlusInj;
2187          l3dHistGenPtVsGenYvsNtracksInj = f3dHistGenPtVsGenYvsNtracksXiPlusInj;
2188          l3dHistGenPtVsGenctauvsCentNat = f3dHistGenPtVsGenctauvsCentXiPlusNat;
2189          l3dHistGenPtVsGenctauvsCentInj = f3dHistGenPtVsGenctauvsCentXiPlusInj;
2190
2191                 // cascades generated within acceptance (cut in pt + theta)
2192          lHistThetaGenCascNat      = fHistThetaGenCascXiPlusNat;
2193          lHistThetaGenCascInj      = fHistThetaGenCascXiPlusInj;
2194          l2dHistGenPtVsGenYFdbl = f2dHistGenPtVsGenYFdblXiPlus;
2195          lHistThetaLambda       = fHistThetaLambdaXiPlus;
2196          lHistThetaBach         = fHistThetaBachXiPlus;
2197          lHistThetaBarDghter    = fHistThetaBarDghterXiPlus;
2198          lHistThetaMesDghter    = fHistThetaMesDghterXiPlus;
2199          lHistPtBach            = fHistPtBachXiPlus;
2200          lHistPtBarDghter       = fHistPtBarDghterXiPlus;
2201          lHistPtMesDghter       = fHistPtMesDghterXiPlus;
2202         break;
2203
2204     case 3: // Omega-
2205          lPdgCodeCasc       =   3334;  //Omega-
2206          lPdgCodeBach       =   -321;  //K-
2207          lPdgCodeLambda     =   3122;  //Lambda0
2208          lPdgCodeDghtMesV0  =   -211;  //Pi-
2209          lPdgCodeDghtBarV0  =   2212;  //Proton
2210
2211                 // any Omega-
2212          lHistEtaGenCasc        = fHistEtaGenCascOmegaMinus;
2213          l3dHistGenPtVsGenYvsCentNat = f3dHistGenPtVsGenYvsCentOmegaMinusNat;
2214          l3dHistGenPtVsGenYvsNtracksNat = f3dHistGenPtVsGenYvsNtracksOmegaMinusNat;
2215          l3dHistGenPtVsGenYvsCentInj = f3dHistGenPtVsGenYvsCentOmegaMinusInj;
2216          l3dHistGenPtVsGenYvsNtracksInj = f3dHistGenPtVsGenYvsNtracksOmegaMinusInj;
2217          l3dHistGenPtVsGenctauvsCentNat = f3dHistGenPtVsGenctauvsCentOmegaMinusNat;
2218          l3dHistGenPtVsGenctauvsCentInj = f3dHistGenPtVsGenctauvsCentOmegaMinusInj;
2219
2220                 // cascades generated within acceptance (cut in pt + theta)
2221          lHistThetaGenCascNat      = fHistThetaGenCascOmegaMinusNat;
2222          lHistThetaGenCascInj      = fHistThetaGenCascOmegaMinusInj;
2223          l2dHistGenPtVsGenYFdbl = f2dHistGenPtVsGenYFdblOmegaMinus;
2224          lHistThetaLambda       = fHistThetaLambdaOmegaMinus;
2225          lHistThetaBach         = fHistThetaBachOmegaMinus;
2226          lHistThetaBarDghter    = fHistThetaBarDghterOmegaMinus;
2227          lHistThetaMesDghter    = fHistThetaMesDghterOmegaMinus;
2228          lHistPtBach            = fHistPtBachOmegaMinus;
2229          lHistPtBarDghter       = fHistPtBarDghterOmegaMinus;
2230          lHistPtMesDghter       = fHistPtMesDghterOmegaMinus;
2231         break;
2232
2233     case 4:  // Omega+
2234          lPdgCodeCasc       =  -3334;  //Omega+
2235          lPdgCodeBach       =    321;  //K+
2236          lPdgCodeLambda     =  -3122;  //AntiLambda0
2237          lPdgCodeDghtMesV0  =    211;  //Pi+
2238          lPdgCodeDghtBarV0  =  -2212;  //AntiProton
2239
2240                 // any Omega+
2241          lHistEtaGenCasc        = fHistEtaGenCascOmegaPlus;
2242          l3dHistGenPtVsGenYvsCentNat = f3dHistGenPtVsGenYvsCentOmegaPlusNat;
2243          l3dHistGenPtVsGenYvsNtracksNat = f3dHistGenPtVsGenYvsNtracksOmegaPlusNat;
2244          l3dHistGenPtVsGenYvsCentInj = f3dHistGenPtVsGenYvsCentOmegaPlusInj;
2245          l3dHistGenPtVsGenYvsNtracksInj = f3dHistGenPtVsGenYvsNtracksOmegaPlusInj;
2246          l3dHistGenPtVsGenctauvsCentNat = f3dHistGenPtVsGenctauvsCentOmegaPlusNat;
2247          l3dHistGenPtVsGenctauvsCentInj = f3dHistGenPtVsGenctauvsCentOmegaPlusInj;
2248
2249
2250                 // cascades generated within acceptance (cut in pt + theta)
2251          lHistThetaGenCascNat      = fHistThetaGenCascOmegaPlusNat;
2252          lHistThetaGenCascInj      = fHistThetaGenCascOmegaPlusInj;
2253          l2dHistGenPtVsGenYFdbl = f2dHistGenPtVsGenYFdblOmegaPlus;
2254          lHistThetaLambda       = fHistThetaLambdaOmegaPlus;
2255          lHistThetaBach         = fHistThetaBachOmegaPlus;
2256          lHistThetaBarDghter    = fHistThetaBarDghterOmegaPlus;
2257          lHistThetaMesDghter    = fHistThetaMesDghterOmegaPlus;
2258          lHistPtBach            = fHistPtBachOmegaPlus;
2259          lHistPtBarDghter       = fHistPtBarDghterOmegaPlus;
2260          lHistPtMesDghter       = fHistPtMesDghterOmegaPlus;
2261         break;
2262
2263   }// end switch cascade
2264
2265
2266   for (Int_t iCurrentLabelStack = 0; iCurrentLabelStack < nNumberOfMCPrimaries; iCurrentLabelStack++) {
2267
2268     Bool_t kIsNaturalPart = kFALSE;
2269     Double_t partEnergy = 0.;
2270     Double_t partPz     = 0.;
2271     Double_t partEta    = 0.;
2272     Double_t partTheta  = 0.;
2273     Double_t partP      = 0.;
2274     Double_t partPt     = 0.;
2275     Double_t partVx     = 0.;
2276     Double_t partVy     = 0.;
2277     Double_t partVz     = 0.;
2278     Double_t bacVx      = 0.;
2279     Double_t bacVy      = 0.;
2280     Double_t bacVz      = 0.;
2281     Double_t partMass   = 0.;
2282
2283     if ( fAnalysisType == "ESD" ) {
2284       TParticle* lCurrentParticle = 0x0;
2285                  lCurrentParticle = lMCstack->Particle( iCurrentLabelStack );
2286       if (!lCurrentParticle) {
2287         Printf("Cascade loop %d - MC TParticle pointer to current stack particle = 0x0 ! Skip ...\n", iCurrentLabelStack );
2288         continue;
2289       }
2290       if (!lMCstack->IsPhysicalPrimary(iCurrentLabelStack)) continue;
2291       if (lCurrentParticle->GetPdgCode() == lPdgCodeCasc) {  // Here !
2292         //cout << "Xi- within loop " << iCurrentLabelStack << "/ " << nNumberOfMCPrimaries << endl;
2293
2294         // -  Xi level ... _____________________________________________________________
2295         TParticle* xiMC = 0x0;
2296                    xiMC = lCurrentParticle;
2297         if (!xiMC) {
2298           Printf("MC TParticle pointer to Cascade = 0x0 ! Skip ...");
2299           continue;
2300         }
2301
2302         kIsNaturalPart = lMCevent->IsFromBGEvent(iCurrentLabelStack);
2303         if (!kIsNaturalPart) {
2304           if (!(lCurrentParticle->GetFirstMother()<0)) kIsNaturalPart = kTRUE; // because there are primaries (ALICE definition) not produced in the collision
2305         }
2306         partEnergy = xiMC->Energy();
2307         partPz = xiMC->Pz();
2308         partEta = xiMC->Eta();
2309         partPt = xiMC->Pt();
2310         partP = xiMC->P();
2311         partTheta = xiMC->Theta();
2312         partMass = xiMC->GetMass();
2313         partVx = xiMC->Vx();
2314         partVy = xiMC->Vy();
2315         partVz = xiMC->Vz();
2316
2317         if (xiMC->GetDaughter(0)>=0) {
2318           TParticle *mcBach = lMCstack->Particle(xiMC->GetDaughter(0));
2319           if (mcBach) {
2320             bacVx  = mcBach->Vx();
2321             bacVy  = mcBach->Vy();
2322             bacVz  = mcBach->Vz();
2323           }
2324         }
2325       } else continue;
2326     } else if ( fAnalysisType == "AOD" ) {
2327
2328       AliAODMCParticle *lCurrentParticleaod = (AliAODMCParticle*) arrayMC->At(iCurrentLabelStack);
2329       if (!lCurrentParticleaod) {
2330         Printf("Cascade loop %d - MC TParticle pointer to current stack particle = 0x0 ! Skip ...\n", iCurrentLabelStack );
2331         continue;
2332       }
2333
2334       if (endOfHijingEvent==-1) {  // For injected MC: determine where HIJING event ends
2335         if ((lCurrentParticleaod->GetStatus()==21)||
2336             ( (lCurrentParticleaod->GetPdgCode() == 443) &&
2337               (lCurrentParticleaod->GetMother() == -1)   &&
2338               (lCurrentParticleaod->GetDaughter(0) ==  (iCurrentLabelStack+1)))) {
2339
2340
2341           endOfHijingEvent = iCurrentLabelStack;
2342
2343         }
2344       }
2345
2346       if (!lCurrentParticleaod->IsPhysicalPrimary()) continue;
2347       if (!(lCurrentParticleaod->PdgCode() == lPdgCodeCasc)) continue;
2348
2349
2350       kIsNaturalPart = kTRUE;
2351       if ((iCurrentLabelStack>=endOfHijingEvent)&&(endOfHijingEvent!=-1)&&(lCurrentParticleaod->GetMother()<0)) kIsNaturalPart = kFALSE;
2352
2353       partEnergy = lCurrentParticleaod->E();
2354       partPz = lCurrentParticleaod->Pz();
2355       partEta = lCurrentParticleaod->Eta();
2356       partP = lCurrentParticleaod->P();
2357       partPt = lCurrentParticleaod->Pt();
2358       partTheta = lCurrentParticleaod->Theta();
2359       partMass = lCurrentParticleaod->M(); // FIXME not sure this works, seems not implemented
2360       partVx = lCurrentParticleaod->Xv();
2361       partVy = lCurrentParticleaod->Yv();
2362       partVz = lCurrentParticleaod->Zv();
2363
2364       if (lCurrentParticleaod->GetDaughter(0)>=0) {
2365         AliAODMCParticle *mcBach = (AliAODMCParticle*) arrayMC->At(lCurrentParticleaod->GetDaughter(0));
2366         if (mcBach) {
2367           bacVx  = mcBach->Xv();
2368           bacVy  = mcBach->Yv();
2369           bacVz  = mcBach->Zv();
2370         }
2371       }
2372
2373     }
2374
2375     ncascperevtot++;
2376         // Fill the first histos : = any generated Xi, not necessarily within the acceptance
2377     Double_t lRapXiMC = 0.5*TMath::Log((partEnergy + partPz) / (partEnergy - partPz +1.e-13));
2378         // Calculate proper time
2379     Double_t lctau = TMath::Sqrt((partVx-bacVx)*(partVx-bacVx)+(partVy-bacVy)*(partVy-bacVy)+(partVz-bacVz)*(partVz-bacVz));
2380     if (partP!=0.)    lctau = lctau*partMass/partP;
2381     else lctau = -1.;
2382
2383     Double_t lRadToDeg = 180.0/TMath::Pi();
2384
2385     // Fill the first histos : = any generated Xi, not necessarily within the acceptance
2386     lHistEtaGenCasc     ->Fill( partEta );
2387
2388     if (kIsNaturalPart) {
2389       l3dHistGenPtVsGenYvsCentNat   ->Fill( partPt, lRapXiMC, lcentrality );
2390       l3dHistGenPtVsGenYvsNtracksNat->Fill( partPt, lRapXiMC, lPrimaryTrackMultiplicity );
2391       if (TMath::Abs(lRapXiMC)<0.5) l3dHistGenPtVsGenctauvsCentNat->Fill( partPt, lctau, lcentrality );
2392       lHistThetaGenCascNat          ->Fill( lRadToDeg * partTheta );
2393     } else {
2394       l3dHistGenPtVsGenYvsCentInj   ->Fill( partPt, lRapXiMC, lcentrality );
2395       l3dHistGenPtVsGenYvsNtracksInj->Fill( partPt, lRapXiMC, lPrimaryTrackMultiplicity );
2396       if (TMath::Abs(lRapXiMC)<0.5) l3dHistGenPtVsGenctauvsCentInj->Fill( partPt, lctau, lcentrality );
2397       lHistThetaGenCascInj          ->Fill( lRadToDeg * partTheta );
2398     }
2399
2400     // Check the emission of particle stays within the acceptance of the detector (cut in theta)
2401     if (fApplyAccCut) { if( partTheta < TMath::Pi()/4.0 || partTheta > 3.0*TMath::Pi()/4.0 ) continue;}
2402
2403     Float_t lambdaTheta = 0.;
2404     Float_t bacTheta    = 0.;
2405     Float_t dghtBarV0Theta = 0.;
2406     Float_t dghtMesV0Theta = 0.;
2407     Float_t bacPt       = 0.;
2408     Float_t dghtBarV0Pt = 0.;
2409     Float_t dghtMesV0Pt = 0.;
2410
2411
2412     if ( fAnalysisType == "ESD" ) {
2413       TParticle* xiMC = lMCstack->Particle( iCurrentLabelStack );
2414       if ( xiMC->GetNDaughters() != 2) continue;
2415       if ( xiMC->GetDaughter(0) < 0 )  continue;
2416       if ( xiMC->GetDaughter(1) < 0 )  continue;
2417
2418       TParticle* lDght0ofXi = lMCstack->Particle(  xiMC->GetDaughter(0) );
2419       TParticle* lDght1ofXi = lMCstack->Particle(  xiMC->GetDaughter(1) );
2420
2421       TParticle* lLambda = 0;
2422       TParticle* lBach   = 0;
2423
2424       // Xi - Case 1
2425       if ( lDght0ofXi->GetPdgCode() == lPdgCodeLambda   &&  // Here !
2426            lDght1ofXi->GetPdgCode() == lPdgCodeBach ){      // Here !
2427
2428         lLambda = lDght0ofXi;
2429         lBach   = lDght1ofXi;
2430       }// end if dghter 0 = Lambda and    dghter 1 = Pi-
2431
2432       // Xi - Case 2
2433         else if ( lDght0ofXi->GetPdgCode() == lPdgCodeBach  &&      // Here !
2434                   lDght1ofXi->GetPdgCode() == lPdgCodeLambda ){     // Here !
2435
2436         lBach   = lDght0ofXi;
2437         lLambda = lDght1ofXi;
2438       }//  end if dghter 0 = Pi-  and   dghter 1 = Lambda
2439
2440         // V0 otherwise - Case 3
2441         else continue;
2442
2443       // Check the emission of particle stays within the acceptance of the detector (cut in pt + theta)
2444       if (fApplyAccCut) {
2445         if ( lLambda->Theta() < TMath::Pi()/4.0  ||    lLambda->Theta() > 3.0*TMath::Pi()/4.0 ) continue;
2446         if( lBach->Theta() < TMath::Pi()/4.0    ||    lBach->Theta() > 3.0*TMath::Pi()/4.0 )   continue;
2447         if( lBach->Pt() < 0.150 ) continue; //FIXME : maybe tuned for Xi but not for K- from Omega ...
2448       }
2449
2450       // -  V0 level ... _____________________________________________________________
2451       TParticle* lDghtBarV0 = 0;
2452       TParticle* lDghtMesV0 = 0;
2453
2454       if( lLambda->GetNDaughters() != 2 )  continue;
2455       if( lLambda->GetDaughter(0) < 0 )    continue;
2456       if( lLambda->GetDaughter(1) < 0 )    continue;
2457
2458
2459       TParticle* lDght0ofLambda = lMCstack->Particle(  lLambda->GetDaughter(0) );
2460       TParticle* lDght1ofLambda = lMCstack->Particle(  lLambda->GetDaughter(1) );
2461
2462       // V0 - Case 1
2463       if ( lDght0ofLambda->GetPdgCode() == lPdgCodeDghtBarV0 &&    // Here !
2464            lDght1ofLambda->GetPdgCode() == lPdgCodeDghtMesV0 ) {    // Here !
2465
2466         lDghtBarV0 = lDght0ofLambda;
2467         lDghtMesV0 = lDght1ofLambda;
2468       }// end if dghter 0 = Proton  and   dghter 1 = Pi-
2469
2470       // V0 - Case 2
2471         else if ( lDght0ofLambda->GetPdgCode() == lPdgCodeDghtMesV0  &&     // Here !
2472                   lDght1ofLambda->GetPdgCode() == lPdgCodeDghtBarV0 ) {      // Here !
2473
2474         lDghtMesV0 = lDght0ofLambda;
2475         lDghtBarV0 = lDght1ofLambda;
2476       }//  end if dghter 0 = Pi-  and   dghter 1 = proton
2477
2478       // V0 otherwise - Case 3
2479         else continue;
2480
2481
2482       // Check the emission of particle stays within the acceptance of the detector
2483       if (fApplyAccCut) {
2484         if( lDghtBarV0->Theta() < TMath::Pi()/4.0  ||  lDghtBarV0->Theta() > 3.0*TMath::Pi()/4.0 ) continue;
2485         if( lDghtMesV0->Theta() < TMath::Pi()/4.0  ||  lDghtMesV0->Theta() > 3.0*TMath::Pi()/4.0 ) continue;
2486
2487         if( lDghtBarV0->Pt() < 0.250 ) continue;
2488         if( lDghtMesV0->Pt() < 0.150 ) continue;
2489       }
2490
2491       lambdaTheta = lLambda->Theta();
2492       bacTheta    = lBach->Theta();
2493       dghtBarV0Theta = lDghtBarV0->Theta();
2494       dghtMesV0Theta = lDghtMesV0->Theta();
2495       bacPt       = lBach->Pt();
2496       dghtBarV0Pt = lDghtBarV0->Pt();
2497       dghtMesV0Pt = lDghtMesV0->Pt();
2498
2499     } else if ( fAnalysisType == "AOD") {
2500
2501       AliAODMCParticle *xiMC = (AliAODMCParticle*) arrayMC->At(iCurrentLabelStack);
2502       if ( xiMC->GetNDaughters() != 2) continue;
2503       if ( xiMC->GetDaughter(0) < 0 )  continue;
2504       if ( xiMC->GetDaughter(1) < 0 )  continue;
2505
2506       AliAODMCParticle* lDght0ofXi = (AliAODMCParticle*) arrayMC->At(  xiMC->GetDaughter(0) );
2507       AliAODMCParticle* lDght1ofXi = (AliAODMCParticle*) arrayMC->At(  xiMC->GetDaughter(1) );
2508
2509       AliAODMCParticle* lLambda = 0;
2510       AliAODMCParticle* lBach   = 0;
2511
2512       // Xi - Case 1
2513       if ( lDght0ofXi->PdgCode() == lPdgCodeLambda   &&  // Here !
2514            lDght1ofXi->PdgCode() == lPdgCodeBach ){      // Here !
2515
2516         lLambda = lDght0ofXi;
2517         lBach   = lDght1ofXi;
2518       }// end if dghter 0 = Lambda and    dghter 1 = Pi-
2519
2520       // Xi - Case 2
2521         else if ( lDght0ofXi->PdgCode() == lPdgCodeBach  &&      // Here !
2522                   lDght1ofXi->PdgCode() == lPdgCodeLambda ){     // Here !
2523
2524         lBach   = lDght0ofXi;
2525         lLambda = lDght1ofXi;
2526       }//  end if dghter 0 = Pi-  and   dghter 1 = Lambda
2527
2528         // V0 otherwise - Case 3
2529         else continue;
2530
2531       // Check the emission of particle stays within the acceptance of the detector (cut in pt + theta)
2532       if (fApplyAccCut) {
2533         if ( lLambda->Theta() < TMath::Pi()/4.0  ||    lLambda->Theta() > 3.0*TMath::Pi()/4.0 ) continue;
2534         if( lBach->Theta() < TMath::Pi()/4.0    ||    lBach->Theta() > 3.0*TMath::Pi()/4.0 )   continue;
2535         if( lBach->Pt() < 0.150 ) continue; //FIXME : maybe tuned for Xi but not for K- from Omega ...
2536       }
2537
2538       // -  V0 level ... _____________________________________________________________
2539       AliAODMCParticle* lDghtBarV0 = 0;
2540       AliAODMCParticle* lDghtMesV0 = 0;
2541
2542       if( lLambda->GetNDaughters() != 2 )  continue;
2543       if( lLambda->GetDaughter(0) < 0 )    continue;
2544       if( lLambda->GetDaughter(1) < 0 )    continue;
2545
2546
2547       AliAODMCParticle* lDght0ofLambda = (AliAODMCParticle*) arrayMC->At(  lLambda->GetDaughter(0) );
2548       AliAODMCParticle* lDght1ofLambda = (AliAODMCParticle*) arrayMC->At(  lLambda->GetDaughter(1) );
2549
2550       // V0 - Case 1
2551       if ( lDght0ofLambda->PdgCode() == lPdgCodeDghtBarV0 &&    // Here !
2552            lDght1ofLambda->PdgCode() == lPdgCodeDghtMesV0 ) {    // Here !
2553
2554         lDghtBarV0 = lDght0ofLambda;
2555         lDghtMesV0 = lDght1ofLambda;
2556       }// end if dghter 0 = Proton  and   dghter 1 = Pi-
2557
2558       // V0 - Case 2
2559         else if ( lDght0ofLambda->PdgCode() == lPdgCodeDghtMesV0  &&     // Here !
2560                   lDght1ofLambda->PdgCode() == lPdgCodeDghtBarV0 ) {      // Here !
2561
2562         lDghtMesV0 = lDght0ofLambda;
2563         lDghtBarV0 = lDght1ofLambda;
2564       }//  end if dghter 0 = Pi-  and   dghter 1 = proton
2565
2566       // V0 otherwise - Case 3
2567         else continue;
2568
2569
2570       // Check the emission of particle stays within the acceptance of the detector
2571       if (fApplyAccCut) {
2572         if( lDghtBarV0->Theta() < TMath::Pi()/4.0  ||  lDghtBarV0->Theta() > 3.0*TMath::Pi()/4.0 ) continue;
2573         if( lDghtMesV0->Theta() < TMath::Pi()/4.0  ||  lDghtMesV0->Theta() > 3.0*TMath::Pi()/4.0 ) continue;
2574
2575         if( lDghtBarV0->Pt() < 0.250 ) continue;
2576         if( lDghtMesV0->Pt() < 0.150 ) continue;
2577       }
2578
2579
2580       lambdaTheta = lLambda->Theta();
2581       bacTheta    = lBach->Theta();
2582       dghtBarV0Theta = lDghtBarV0->Theta();
2583       dghtMesV0Theta = lDghtMesV0->Theta();
2584       bacPt       = lBach->Pt();
2585       dghtBarV0Pt = lDghtBarV0->Pt();
2586       dghtMesV0Pt = lDghtMesV0->Pt();
2587
2588     }
2589
2590     // - Just to know which file is currently open : locate the file containing Xi
2591                 //cout << "Name of the file containing generated Xi :" <<  CurrentFileName()
2592                 //
2593     // - Filling histos for findable cascades... _________________________________________________________________
2594
2595     // - Fill theta histos for Lambda and Bach
2596     lHistThetaLambda        ->Fill( lRadToDeg * lambdaTheta );
2597     lHistThetaBach          ->Fill( lRadToDeg * bacTheta );
2598
2599     // - Fill theta histos for V0 daughters
2600     lHistThetaBarDghter     ->Fill( lRadToDeg * dghtBarV0Theta );
2601     lHistThetaMesDghter     ->Fill( lRadToDeg * dghtMesV0Theta );
2602
2603     // - Fill pt histos.
2604     lHistPtBach             ->Fill( bacPt );
2605     lHistPtBarDghter        ->Fill( dghtBarV0Pt );
2606     lHistPtMesDghter        ->Fill( dghtMesV0Pt );
2607     //  if(iCascType == 1) Printf("Xi- current index = %n ", iCurrentLabelStack);
2608
2609     l2dHistGenPtVsGenYFdbl  ->Fill( partPt, lRapXiMC );
2610
2611     ncascperev++;
2612
2613
2614   }// This is the end of the loop on primaries
2615
2616   if (iCascType == 1) {
2617     fHistnXiMinusPerEv->Fill(ncascperev);
2618     fHistnXiMinusPerEvTot->Fill(ncascperevtot);
2619 //       Printf("N xi-tot = %n N xi-acc = %n\n", ncascperevtot, ncascperev);
2620   }
2621   if (iCascType == 2) {
2622     fHistnXiPlusPerEv->Fill(ncascperev);
2623     fHistnXiPlusPerEvTot->Fill(ncascperevtot);
2624   }
2625   if (iCascType == 3) {
2626     fHistnOmegaMinusPerEv->Fill(ncascperev);
2627     fHistnOmegaMinusPerEvTot->Fill(ncascperevtot);
2628   }
2629   if (iCascType == 4) {
2630     fHistnOmegaPlusPerEv->Fill(ncascperev);
2631     fHistnOmegaPlusPerEvTot->Fill(ncascperevtot);
2632   }
2633
2634
2635 // - Re-initialisation of the local TH1F pointers
2636 lHistEtaGenCasc         = 0x0;
2637 l3dHistGenPtVsGenYvsCentNat = 0x0;
2638 l3dHistGenPtVsGenYvsNtracksNat = 0x0;
2639 l3dHistGenPtVsGenYvsCentInj = 0x0;
2640 l3dHistGenPtVsGenYvsNtracksInj = 0x0;
2641 l3dHistGenPtVsGenctauvsCentNat = 0x0;
2642 l3dHistGenPtVsGenctauvsCentInj = 0x0;
2643
2644 lHistThetaGenCascNat    = 0x0;
2645 lHistThetaGenCascInj    = 0x0;
2646 l2dHistGenPtVsGenYFdbl  = 0x0;
2647 lHistThetaLambda        = 0x0;
2648 lHistThetaBach          = 0x0;
2649 lHistThetaBarDghter     = 0x0;
2650 lHistThetaMesDghter     = 0x0;
2651 lHistPtBach             = 0x0;
2652 lHistPtBarDghter        = 0x0;
2653 lHistPtMesDghter        = 0x0;
2654
2655 } // end of loop over the different types of cascades (Xi-+, Omega-+)
2656
2657
2658
2659 //__________________________________________________________________________
2660 // Part 3 - Loop over the reconstructed candidates
2661
2662 Int_t nAssoXiMinus = 0;
2663 Int_t nAssoXiPlus = 0;
2664 Int_t nAssoOmegaMinus = 0;
2665 Int_t nAssoOmegaPlus = 0;
2666
2667
2668 Int_t lPosTPCClusters   = 0;
2669 Int_t lNegTPCClusters   = 0;
2670 Int_t lBachTPCClusters  = 0;
2671
2672
2673 // Double_t lChi2Xi         = -1. ;
2674 Double_t lDcaXiDaughters            = -1. ;
2675 Double_t lDcaBachToPrimVertexXi     = -1. ;
2676 Double_t lXiCosineOfPointingAngle   = -1. ;
2677 Double_t lPosXi[3]                  = { -1000.0, -1000.0, -1000.0 };
2678 Double_t lXiRadius                  = -1000. ;
2679
2680 Double_t lInvMassLambdaAsCascDghter = 0.;
2681 Double_t lDcaV0DaughtersXi          = -1.;
2682 // Double_t lV0Chi2Xi               = -1. ;
2683 Double_t lV0CosineOfPointingAngle = -1.;
2684 Double_t lV0toXiCosineOfPointingAngle   = -1.;
2685 Double_t lPosV0Xi[3]                = { -1000. , -1000., -1000. }; // Position of VO coming from cascade
2686 Double_t lV0RadiusXi                = -1000.;
2687 Double_t lDcaV0ToPrimVertexXi       = -1.;
2688 Double_t lDcaPosToPrimVertexXi      = -1.;
2689 Double_t lDcaNegToPrimVertexXi      = -1.;
2690
2691 Double_t lChargeXi                  = -1.;
2692
2693 Double_t lV0mom = -1000.;
2694
2695 Double_t lmcPt = -1.;
2696 Double_t lmcRapCasc = -1.;
2697 Double_t lmcEta = -1000.;
2698 Double_t lmcTransvRadius = -1000.;
2699
2700
2701 Double_t lrecoPt = -100.;
2702 Double_t lrecoTransvRadius = -1000.;
2703
2704 Double_t lDeltaPhiMcReco = -1.;
2705
2706
2707 Double_t lBachTransvMom = 0.;
2708 Double_t lpTrackTransvMom  = 0.;
2709 Double_t lnTrackTransvMom  = 0.;
2710
2711 Double_t lmcPtPosV0Dghter = -100.;
2712 Double_t lmcPtNegV0Dghter = -100.;
2713 Double_t lrecoP = -100.;
2714
2715 Double_t   lmcPtBach = -100.;
2716
2717 Double_t cascadeMass = 0.;
2718
2719
2720
2721 for (Int_t iXi = 0; iXi < ncascades; iXi++) {// This is the begining of the Cascade loop
2722
2723
2724       Bool_t   lIsPosInXiProton      = kFALSE;
2725       Bool_t   lIsPosInXiPion        = kFALSE;
2726       Bool_t   lIsPosInOmegaProton   = kFALSE;
2727       Bool_t   lIsPosInOmegaPion     = kFALSE;
2728
2729       Bool_t   lIsNegInXiProton      = kFALSE;
2730       Bool_t   lIsNegInXiPion        = kFALSE;
2731       Bool_t   lIsNegInOmegaProton   = kFALSE;
2732       Bool_t   lIsNegInOmegaPion     = kFALSE;
2733
2734       Bool_t   lIsBachelorKaon       = kFALSE;
2735       Bool_t   lIsBachelorPion       = kFALSE;
2736
2737       Bool_t   lIsBachelorKaonForTPC = kFALSE;
2738       Bool_t   lIsBachelorPionForTPC = kFALSE;
2739       Bool_t   lIsNegPionForTPC      = kFALSE;
2740       Bool_t   lIsPosPionForTPC      = kFALSE;
2741       Bool_t   lIsNegProtonForTPC    = kFALSE;
2742       Bool_t   lIsPosProtonForTPC    = kFALSE;
2743
2744       // Combined PID
2745       // Reasonable guess for the priors for the cascade track sample (e-, mu, pi, K, p)
2746       Double_t lPriorsGuessXi[5]    = {0, 0, 2, 0, 1};
2747       Double_t lPriorsGuessOmega[5] = {0, 0, 1, 1, 1};
2748
2749       Double_t ppionBach = 0.0, pkaonBach = 0.0;
2750
2751       Bool_t   lIsBachelorMCPiMinus  = kFALSE;
2752       Bool_t   lIsBachelorMCPiPlus   = kFALSE;
2753       Bool_t   lIsBachelorMCKMinus   = kFALSE;
2754       Bool_t   lIsBachelorMCKPlus    = kFALSE;
2755
2756       Double_t lInvMassXiMinus    = 0.;
2757       Double_t lInvMassXiPlus     = 0.;
2758       Double_t lInvMassOmegaMinus = 0.;
2759       Double_t lInvMassOmegaPlus  = 0.;
2760
2761       Bool_t lAssoXiMinus    = kFALSE;
2762       Bool_t lAssoXiPlus     = kFALSE;
2763       Bool_t lAssoOmegaMinus = kFALSE;
2764       Bool_t lAssoOmegaPlus  = kFALSE;
2765
2766       Bool_t kIsNaturalPart = kFALSE;
2767
2768       Float_t  etaBach = 0.;
2769       Float_t  etaPos  = 0.;
2770       Float_t  etaNeg  = 0.;
2771
2772
2773       if ( fAnalysisType == "ESD" ) {
2774         AliESDcascade *xiESD = lESDevent->GetCascade(iXi);
2775         if (!xiESD) continue;
2776
2777         // - Step II.1 : Connection to daughter tracks of the current cascade
2778         //-------------
2779
2780                 UInt_t lIdxPosXi        = (UInt_t) TMath::Abs( xiESD->GetPindex() );
2781                 UInt_t lIdxNegXi        = (UInt_t) TMath::Abs( xiESD->GetNindex() );
2782                 UInt_t lBachIdx         = (UInt_t) TMath::Abs( xiESD->GetBindex() );
2783                 // abs value not needed ; the index should always be positive (!= label ...)
2784
2785
2786         // FIXME : rejection of a double use of a daughter track (nothing but just a crosscheck of what is done in the cascade vertexer)
2787         if(lBachIdx == lIdxNegXi) {
2788                 AliWarning("Pb / Idx(Bach. track) = Idx(Neg. track) ... continue!"); continue;
2789         }
2790         if(lBachIdx == lIdxPosXi) {
2791                 AliWarning("Pb / Idx(Bach. track) = Idx(Pos. track) ... continue!"); continue;
2792         }
2793
2794         AliESDtrack *pTrackXi           = lESDevent->GetTrack( lIdxPosXi );
2795         AliESDtrack *nTrackXi           = lESDevent->GetTrack( lIdxNegXi );
2796         AliESDtrack *bachTrackXi        = lESDevent->GetTrack( lBachIdx  );
2797         if (!pTrackXi || !nTrackXi || !bachTrackXi ) {
2798                 Printf("ERROR: Could not retrieve one of the 3 daughter tracks of the cascade ...");
2799                 continue;
2800         }
2801
2802         lPosTPCClusters   = pTrackXi->GetTPCNcls();
2803         lNegTPCClusters   = nTrackXi->GetTPCNcls();
2804         lBachTPCClusters  = bachTrackXi->GetTPCNcls();
2805
2806                 // FIXME : rejection of a poor quality tracks
2807         if(fkQualityCutTPCrefit){
2808                 // 1 - Poor quality related to TPCrefit
2809                 ULong_t pStatus    = pTrackXi->GetStatus();
2810                 ULong_t nStatus    = nTrackXi->GetStatus();
2811                 ULong_t bachStatus = bachTrackXi->GetStatus();
2812                 if ((pStatus&AliESDtrack::kTPCrefit)    == 0) { AliWarning("Pb / V0 Pos. track has no TPCrefit ... continue!"); continue; }
2813                 if ((nStatus&AliESDtrack::kTPCrefit)    == 0) { AliWarning("Pb / V0 Neg. track has no TPCrefit ... continue!"); continue; }
2814                 if ((bachStatus&AliESDtrack::kTPCrefit) == 0) { AliWarning("Pb / Bach.   track has no TPCrefit ... continue!"); continue; }
2815         }
2816         if(fkQualityCutnTPCcls){
2817                 // 2 - Poor quality related to TPC clusters
2818                 if(lPosTPCClusters  < fMinnTPCcls) { AliWarning("Pb / V0 Pos. track has less than n TPC clusters ... continue!"); continue; }
2819                 if(lNegTPCClusters  < fMinnTPCcls) { AliWarning("Pb / V0 Neg. track has less than n TPC clusters ... continue!"); continue; }
2820                 if(lBachTPCClusters < fMinnTPCcls) { AliWarning("Pb / Bach.   track has less than n TPC clusters ... continue!"); continue; }
2821         }
2822
2823         etaPos  = pTrackXi->Eta();
2824         etaNeg  = nTrackXi->Eta();
2825         etaBach = bachTrackXi->Eta();
2826
2827
2828         // - Step II.2 : Info over reconstructed cascades
2829         //-------------
2830
2831
2832         Double_t lV0quality = 0.;
2833
2834         if( bachTrackXi->Charge() < 0 ) {
2835                 lV0quality = 0.;
2836                 xiESD->ChangeMassHypothesis(lV0quality , 3312);
2837                         // Calculate the effective mass of the Xi- candidate.
2838                         // pdg code 3312 = Xi-
2839                 lInvMassXiMinus = xiESD->GetEffMassXi();
2840
2841                 lV0quality = 0.;
2842                 xiESD->ChangeMassHypothesis(lV0quality , 3334);
2843                         // Calculate the effective mass of the Xi- candidate.
2844                         // pdg code 3334 = Omega-
2845                 lInvMassOmegaMinus = xiESD->GetEffMassXi();
2846
2847                 lV0quality = 0.;
2848                 xiESD->ChangeMassHypothesis(lV0quality , 3312);         // Back to default hyp.
2849
2850         }
2851
2852         if( bachTrackXi->Charge() >  0 ){
2853                 lV0quality = 0.;
2854                 xiESD->ChangeMassHypothesis(lV0quality , -3312);
2855                         // Calculate the effective mass of the Xi+ candidate.
2856                         // pdg code -3312 = Xi+
2857                 lInvMassXiPlus = xiESD->GetEffMassXi();
2858
2859                 lV0quality = 0.;
2860                 xiESD->ChangeMassHypothesis(lV0quality , -3334);
2861                         // Calculate the effective mass of the Xi+ candidate.
2862                         // pdg code -3334  = Omega+
2863                 lInvMassOmegaPlus = xiESD->GetEffMassXi();
2864
2865                 lV0quality = 0.;
2866                 xiESD->ChangeMassHypothesis(lV0quality , -3312);        // Back to "default" hyp.
2867         }
2868
2869
2870         //lChi2Xi                       = xiESD->GetChi2Xi();
2871         lDcaXiDaughters                 = xiESD->GetDcaXiDaughters();
2872         lDcaBachToPrimVertexXi          = TMath::Abs( bachTrackXi->GetD( lBestPrimaryVtxPos[0],
2873                                                                          lBestPrimaryVtxPos[1],
2874                                                                          lMagneticField  ) );
2875                                         // NOTE : AliExternalTrackParam::GetD returns an algebraic value
2876         lXiCosineOfPointingAngle        = xiESD->GetCascadeCosineOfPointingAngle( lBestPrimaryVtxPos[0],
2877                                                                                   lBestPrimaryVtxPos[1],
2878                                                                                   lBestPrimaryVtxPos[2] );
2879                                         // Take care : the best available vertex should be used (like in AliCascadeVertexer)
2880         xiESD->GetXYZcascade( lPosXi[0],  lPosXi[1], lPosXi[2] );
2881         lInvMassLambdaAsCascDghter      = xiESD->GetEffMass();
2882                                         // This value shouldn't change, whatever the working hyp. is : Xi-, Xi+, Omega-, Omega+
2883         lDcaV0DaughtersXi               = xiESD->GetDcaV0Daughters();
2884         // lV0Chi2Xi                    = xiESD->GetChi2V0();
2885         lV0toXiCosineOfPointingAngle    = xiESD->GetV0CosineOfPointingAngle( lPosXi[0],
2886                                                                              lPosXi[1],
2887                                                                              lPosXi[2] );
2888
2889         lV0CosineOfPointingAngle      = xiESD->GetV0CosineOfPointingAngle( lBestPrimaryVtxPos[0],
2890                                                                            lBestPrimaryVtxPos[1],
2891                                                                            lBestPrimaryVtxPos[2]);
2892
2893         //if (lV0toXiCosineOfPointingAngle==1.) cout << "Cosine V0 PA wrt Xi pos ==1!" <<endl;
2894         //else if (lV0toXiCosineOfPointingAngle>1.) cout <<"Cosine V0 PA wrt Xi pos >1!"<<endl;
2895         xiESD->GetXYZ( lPosV0Xi[0],  lPosV0Xi[1], lPosV0Xi[2] );
2896
2897         lDcaV0ToPrimVertexXi            = xiESD->GetD( lBestPrimaryVtxPos[0],
2898                                                 lBestPrimaryVtxPos[1],
2899                                                 lBestPrimaryVtxPos[2] );
2900
2901         lDcaPosToPrimVertexXi           = TMath::Abs( pTrackXi   ->GetD( lBestPrimaryVtxPos[0],
2902                                                                          lBestPrimaryVtxPos[1],
2903                                                                          lMagneticField  )     );
2904
2905         lDcaNegToPrimVertexXi           = TMath::Abs( nTrackXi   ->GetD( lBestPrimaryVtxPos[0],
2906                                                                          lBestPrimaryVtxPos[1],
2907                                                                          lMagneticField  )     );
2908
2909         lChargeXi = xiESD->Charge();
2910
2911         // - Step II.3 : PID info
2912
2913         //-------------
2914
2915
2916         // 3.1 - PID Information
2917
2918
2919         // Combined VO-positive-daughter PID
2920         AliPID pPidXi;          pPidXi.SetPriors(    lPriorsGuessXi    );
2921         AliPID pPidOmega;       pPidOmega.SetPriors( lPriorsGuessOmega );
2922
2923         if( pTrackXi->IsOn(AliESDtrack::kESDpid) ){  // Combined PID exists
2924                 Double_t r[10] = {0.}; pTrackXi->GetESDpid(r);
2925                 pPidXi.SetProbabilities(r);
2926                 pPidOmega.SetProbabilities(r);
2927
2928                 // Check if the V0 positive track is a proton (case for Xi-)
2929                 Double_t pproton = pPidXi.GetProbability(AliPID::kProton);
2930                 if (pproton > pPidXi.GetProbability(AliPID::kElectron) &&
2931                     pproton > pPidXi.GetProbability(AliPID::kMuon)     &&
2932                     pproton > pPidXi.GetProbability(AliPID::kPion)     &&
2933                     pproton > pPidXi.GetProbability(AliPID::kKaon)     )     lIsPosInXiProton = kTRUE;
2934
2935                 // Check if the V0 positive track is a pi+ (case for Xi+)
2936                 Double_t ppion = pPidXi.GetProbability(AliPID::kPion);
2937                 if (ppion > pPidXi.GetProbability(AliPID::kElectron) &&
2938                     ppion > pPidXi.GetProbability(AliPID::kMuon)     &&
2939                     ppion > pPidXi.GetProbability(AliPID::kKaon)     &&
2940                     ppion > pPidXi.GetProbability(AliPID::kProton)   )     lIsPosInXiPion = kTRUE;
2941
2942
2943                 // Check if the V0 positive track is a proton (case for Omega-)
2944                 pproton = 0.;
2945                     pproton = pPidOmega.GetProbability(AliPID::kProton);
2946                 if (pproton > pPidOmega.GetProbability(AliPID::kElectron) &&
2947                     pproton > pPidOmega.GetProbability(AliPID::kMuon)     &&
2948                     pproton > pPidOmega.GetProbability(AliPID::kPion)     &&
2949                     pproton > pPidOmega.GetProbability(AliPID::kKaon)     )  lIsPosInOmegaProton = kTRUE;
2950
2951                 // Check if the V0 positive track is a pi+ (case for Omega+)
2952                 ppion = 0.;
2953                     ppion = pPidOmega.GetProbability(AliPID::kPion);
2954                 if (ppion > pPidOmega.GetProbability(AliPID::kElectron) &&
2955                     ppion > pPidOmega.GetProbability(AliPID::kMuon)     &&
2956                     ppion > pPidOmega.GetProbability(AliPID::kKaon)     &&
2957                     ppion > pPidOmega.GetProbability(AliPID::kProton)   )    lIsPosInOmegaPion = kTRUE;
2958
2959         }// end if V0 positive track with existing combined PID
2960
2961
2962         // Combined VO-negative-daughter PID
2963         AliPID nPidXi;          nPidXi.SetPriors(    lPriorsGuessXi    );
2964         AliPID nPidOmega;       nPidOmega.SetPriors( lPriorsGuessOmega );
2965
2966         if( nTrackXi->IsOn(AliESDtrack::kESDpid) ) {  // Combined PID exists
2967                 Double_t r[10] = {0.}; nTrackXi->GetESDpid(r);
2968                 nPidXi.SetProbabilities(r);
2969                 nPidOmega.SetProbabilities(r);
2970
2971                 // Check if the V0 negative track is a pi- (case for Xi-)
2972                 Double_t ppion = nPidXi.GetProbability(AliPID::kPion);
2973                 if (ppion > nPidXi.GetProbability(AliPID::kElectron) &&
2974                     ppion > nPidXi.GetProbability(AliPID::kMuon)     &&
2975                     ppion > nPidXi.GetProbability(AliPID::kKaon)     &&
2976                     ppion > nPidXi.GetProbability(AliPID::kProton)   )     lIsNegInXiPion = kTRUE;
2977
2978                 // Check if the V0 negative track is an anti-proton (case for Xi+)
2979                 Double_t pproton = nPidXi.GetProbability(AliPID::kProton);
2980                 if (pproton > nPidXi.GetProbability(AliPID::kElectron) &&
2981                     pproton > nPidXi.GetProbability(AliPID::kMuon)     &&
2982                     pproton > nPidXi.GetProbability(AliPID::kPion)     &&
2983                     pproton > nPidXi.GetProbability(AliPID::kKaon)     )     lIsNegInXiProton = kTRUE;
2984
2985                 // Check if the V0 negative track is a pi- (case for Omega-)
2986                 ppion = 0.;
2987                     ppion = nPidOmega.GetProbability(AliPID::kPion);
2988                 if (ppion > nPidOmega.GetProbability(AliPID::kElectron) &&
2989                     ppion > nPidOmega.GetProbability(AliPID::kMuon)     &&
2990                     ppion > nPidOmega.GetProbability(AliPID::kKaon)     &&
2991                     ppion > nPidOmega.GetProbability(AliPID::kProton)   )    lIsNegInOmegaPion = kTRUE;
2992
2993                 // Check if the V0 negative track is an anti-proton (case for Omega+)
2994                 pproton = 0.;
2995                     pproton = nPidOmega.GetProbability(AliPID::kProton);
2996                 if (pproton > nPidOmega.GetProbability(AliPID::kElectron) &&
2997                     pproton > nPidOmega.GetProbability(AliPID::kMuon)     &&
2998                     pproton > nPidOmega.GetProbability(AliPID::kPion)     &&
2999                     pproton > nPidOmega.GetProbability(AliPID::kKaon)     )  lIsNegInOmegaProton = kTRUE;
3000
3001         }// end if V0 negative track with existing combined PID
3002
3003
3004         // Combined bachelor PID
3005         AliPID bachPidXi;       bachPidXi.SetPriors(    lPriorsGuessXi    );
3006         AliPID bachPidOmega;    bachPidOmega.SetPriors( lPriorsGuessOmega );
3007
3008
3009         if ( bachTrackXi->IsOn(AliESDtrack::kESDpid) ) {  // Combined PID exists
3010                 Double_t r[10] = {0.}; bachTrackXi->GetESDpid(r);
3011                 bachPidXi.SetProbabilities(r);
3012                 bachPidOmega.SetProbabilities(r);
3013                 // Check if the bachelor track is a pion
3014                     ppionBach = bachPidXi.GetProbability(AliPID::kPion);
3015                 if (ppionBach > bachPidXi.GetProbability(AliPID::kElectron) &&
3016                     ppionBach > bachPidXi.GetProbability(AliPID::kMuon)     &&
3017                     ppionBach > bachPidXi.GetProbability(AliPID::kKaon)     &&
3018                     ppionBach > bachPidXi.GetProbability(AliPID::kProton)   )     lIsBachelorPion = kTRUE;
3019                 // Check if the bachelor track is a kaon
3020                     pkaonBach = bachPidOmega.GetProbability(AliPID::kKaon);
3021                 if (pkaonBach > bachPidOmega.GetProbability(AliPID::kElectron) &&
3022                     pkaonBach > bachPidOmega.GetProbability(AliPID::kMuon)     &&
3023                     pkaonBach > bachPidOmega.GetProbability(AliPID::kPion)     &&
3024                     pkaonBach > bachPidOmega.GetProbability(AliPID::kProton)   )  lIsBachelorKaon = kTRUE;
3025         }// end if bachelor track with existing combined PID
3026
3027
3028         // 3.1.B - TPC PID : 4-sigma bands on Bethe-Bloch curve
3029
3030         // Bachelor
3031         if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( bachTrackXi,AliPID::kKaon)) < 4) lIsBachelorKaonForTPC = kTRUE;
3032         if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( bachTrackXi,AliPID::kPion)) < 4) lIsBachelorPionForTPC = kTRUE;
3033
3034         // Negative V0 daughter
3035         if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( nTrackXi,AliPID::kPion   )) < 4) lIsNegPionForTPC   = kTRUE;
3036         if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( nTrackXi,AliPID::kProton )) < 4) lIsNegProtonForTPC = kTRUE;
3037
3038         // Positive V0 daughter
3039         if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( pTrackXi,AliPID::kPion   )) < 4) lIsPosPionForTPC   = kTRUE;
3040         if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( pTrackXi,AliPID::kProton )) < 4) lIsPosProtonForTPC = kTRUE;
3041
3042 /*
3043         const AliExternalTrackParam *pInnerWallTrackXi    = pTrackXi    ->GetInnerParam(); // Do not use GetTPCInnerWall
3044         const AliExternalTrackParam *nInnerWallTrackXi    = nTrackXi    ->GetInnerParam();
3045         const AliExternalTrackParam *bachInnerWallTrackXi = bachTrackXi ->GetInnerParam();
3046         if(pInnerWallTrackXi && nInnerWallTrackXi && bachInnerWallTrackXi ){
3047
3048                 Double_t pMomInnerWall    = pInnerWallTrackXi   ->GetP();
3049                 Double_t nMomInnerWall    = nInnerWallTrackXi   ->GetP();
3050                 Double_t bachMomInnerWall = bachInnerWallTrackXi->GetP();
3051
3052                 // Bachelor
3053                 if (TMath::Abs(fESDpid->NumberOfSigmasTPC( bachTrackXi,AliPID::kPion)) < 3)                              lIsBachelorPionForTPC = kTRUE;
3054                 if (bachMomInnerWall < 0.350  && TMath::Abs(fESDpid->NumberOfSigmasTPC( bachTrackXi,AliPID::kKaon)) < 5) lIsBachelorKaonForTPC = kTRUE;
3055                 if (bachMomInnerWall > 0.350  && TMath::Abs(fESDpid->NumberOfSigmasTPC( bachTrackXi,AliPID::kKaon)) < 3) lIsBachelorKaonForTPC = kTRUE;
3056
3057                 // Negative V0 daughter
3058                 if (TMath::Abs(fESDpid->NumberOfSigmasTPC( nTrackXi,AliPID::kPion   )) < 3  )                           lIsNegPionForTPC   = kTRUE;
3059                 if (nMomInnerWall < 0.6  && TMath::Abs(fESDpid->NumberOfSigmasTPC( nTrackXi,AliPID::kProton ) ) < 5 )   lIsNegProtonForTPC = kTRUE;
3060                 if (nMomInnerWall > 0.6  && TMath::Abs(fESDpid->NumberOfSigmasTPC( nTrackXi,AliPID::kProton ) ) < 3 )   lIsNegProtonForTPC = kTRUE;
3061
3062                 // Positive V0 daughter
3063                 if (TMath::Abs(fESDpid->NumberOfSigmasTPC( pTrackXi,AliPID::kPion   )) < 3 )                            lIsPosPionForTPC   = kTRUE;
3064                 if (pMomInnerWall < 0.6  && TMath::Abs(fESDpid->NumberOfSigmasTPC( pTrackXi,AliPID::kProton )) < 5)     lIsPosProtonForTPC = kTRUE;
3065                 if (pMomInnerWall > 0.6  && TMath::Abs(fESDpid->NumberOfSigmasTPC( pTrackXi,AliPID::kProton )) < 3)     lIsPosProtonForTPC = kTRUE;
3066
3067         }
3068 */
3069
3070         // 3.2 - PID proba Vs Pt(Bach)
3071         Int_t      lblBachForPID  = (Int_t) TMath::Abs( bachTrackXi->GetLabel() );
3072         TParticle* mcBachForPID   = lMCstack->Particle( lblBachForPID );
3073         lmcPtBach      = mcBachForPID->Pt();
3074
3075         // 3.3 - MC perfect PID
3076
3077         if( mcBachForPID->GetPdgCode() == -211) lIsBachelorMCPiMinus = kTRUE;
3078         if( mcBachForPID->GetPdgCode() ==  211) lIsBachelorMCPiPlus  = kTRUE;
3079         if( mcBachForPID->GetPdgCode() == -321) lIsBachelorMCKMinus  = kTRUE;
3080         if( mcBachForPID->GetPdgCode() ==  321) lIsBachelorMCKPlus   = kTRUE;
3081
3082
3083
3084         // - Step II.4 : MC association (care : lots of "continue;" below this line)
3085         //-------------
3086
3087         if(fDebug > 5)
3088                 cout    << "MC EventNumber : " << lMCevent->Header()->GetEvent()
3089                         << " / MC event Number in Run : " << lMCevent->Header()->GetEventNrInRun() << endl;
3090
3091         // - Step 4.1 : level of the V0 daughters
3092
3093         Int_t lblPosV0Dghter = (Int_t) TMath::Abs( pTrackXi->GetLabel() );
3094                 // Abs value = needed ! question of quality track association ... (negative label when at least one cluster in the track is from a different particle)
3095         Int_t lblNegV0Dghter = (Int_t) TMath::Abs( nTrackXi->GetLabel() );
3096         TParticle* mcPosV0Dghter = lMCstack->Particle( lblPosV0Dghter );
3097         TParticle* mcNegV0Dghter = lMCstack->Particle( lblNegV0Dghter );
3098
3099
3100         // - Step 4.2 : level of the Xi daughters
3101
3102         Int_t lblMotherPosV0Dghter = mcPosV0Dghter->GetFirstMother() ;
3103         Int_t lblMotherNegV0Dghter = mcNegV0Dghter->GetFirstMother();
3104
3105                 if( lblMotherPosV0Dghter != lblMotherNegV0Dghter) continue; // same mother
3106                 if( lblMotherPosV0Dghter < 0 ) continue; // this particle is primary, no mother
3107                 if( lblMotherNegV0Dghter < 0 ) continue;
3108
3109
3110                 // mothers = Lambda candidate ... a priori
3111
3112         TParticle* mcMotherPosV0Dghter = lMCstack->Particle( lblMotherPosV0Dghter );
3113         TParticle* mcMotherNegV0Dghter = lMCstack->Particle( lblMotherNegV0Dghter );  // MN: redundant?? already checked that labels are the same...-->same part from stack
3114
3115         Int_t      lblBach  = (Int_t) TMath::Abs( bachTrackXi->GetLabel() );
3116         TParticle* mcBach   = lMCstack->Particle( lblBach );
3117
3118
3119         // - Step 4.3 : level of Xi candidate
3120
3121         Int_t lblGdMotherPosV0Dghter =   mcMotherPosV0Dghter->GetFirstMother() ;
3122         Int_t lblGdMotherNegV0Dghter =   mcMotherNegV0Dghter->GetFirstMother() ;
3123
3124                 if( lblGdMotherPosV0Dghter != lblGdMotherNegV0Dghter ) continue;
3125                 if( lblGdMotherPosV0Dghter < 0 ) continue; // primary lambda ...
3126                 if( lblGdMotherNegV0Dghter < 0 ) continue; // primary lambda ...
3127
3128                 // Gd mothers = Xi candidate ... a priori
3129
3130         TParticle* mcGdMotherPosV0Dghter = lMCstack->Particle( lblGdMotherPosV0Dghter );
3131         TParticle* mcGdMotherNegV0Dghter = lMCstack->Particle( lblGdMotherNegV0Dghter );
3132
3133         Int_t lblMotherBach = (Int_t) TMath::Abs( mcBach->GetFirstMother()  );
3134
3135         if( lblMotherBach != lblGdMotherPosV0Dghter ) continue; //same mother for bach and V0 daughters
3136
3137         TParticle* mcMotherBach = lMCstack->Particle( lblMotherBach );
3138
3139         // Check if cascade is primary
3140
3141         if (!(lMCstack->IsPhysicalPrimary(lblMotherBach))) continue;
3142
3143         kIsNaturalPart = lMCevent->IsFromBGEvent(lblMotherBach); // if kTRUE the particle is a natural one
3144         if (!kIsNaturalPart) {
3145           if (!(mcMotherBach->GetFirstMother()<0)) kIsNaturalPart = kTRUE;
3146         }
3147
3148         // - Step 4.4 : Manage boolean for association
3149
3150         if( mcMotherBach                ->GetPdgCode() ==   3312 &&
3151             mcGdMotherPosV0Dghter       ->GetPdgCode() ==   3312 &&
3152             mcGdMotherNegV0Dghter       ->GetPdgCode() ==   3312)    {  lAssoXiMinus = kTRUE;
3153                                                                         cascadeMass = 1.321;
3154                                                                       //  Printf("Xi- asso current index = %n ", lblGdMotherPosV0Dghter);
3155                                                                         nAssoXiMinus++; }
3156
3157         else if( mcMotherBach           ->GetPdgCode() ==  -3312 &&
3158             mcGdMotherPosV0Dghter       ->GetPdgCode() ==  -3312 &&
3159             mcGdMotherNegV0Dghter       ->GetPdgCode() ==  -3312)    {  lAssoXiPlus = kTRUE;
3160                                                                         cascadeMass = 1.321;
3161                                                                         nAssoXiPlus++; }
3162
3163         else if( mcMotherBach           ->GetPdgCode() ==   3334 &&
3164             mcGdMotherPosV0Dghter       ->GetPdgCode() ==   3334 &&
3165             mcGdMotherNegV0Dghter       ->GetPdgCode() ==   3334)    {  lAssoOmegaMinus = kTRUE;
3166                                                                         cascadeMass = 1.672;
3167                                                                         nAssoOmegaMinus++; }
3168
3169         else if( mcMotherBach           ->GetPdgCode() ==  -3334 &&
3170             mcGdMotherPosV0Dghter       ->GetPdgCode() ==  -3334 &&
3171             mcGdMotherNegV0Dghter       ->GetPdgCode() ==  -3334)    {  lAssoOmegaPlus = kTRUE;
3172                                                                         cascadeMass = 1.672;
3173                                                                         nAssoOmegaPlus++; }
3174
3175
3176
3177         // If a proper association  exists ...
3178
3179         if(fDebug > 4){
3180                 cout << "XiMinus    = " << lAssoXiMinus    << endl;
3181                 cout << "XiPlus     = " << lAssoXiPlus     << endl;
3182                 cout << "OmegaMinus = " << lAssoOmegaMinus << endl;
3183                 cout << "OmegaPlus  = " << lAssoOmegaPlus  << endl
3184                      << "----"          << endl;
3185         }
3186
3187
3188         if(fDebug > 5){
3189                 cout << endl;
3190                 cout << "- V0 daughters - " << endl;
3191                 cout << "     + V0 Pos. / Label : " << lblPosV0Dghter
3192                 << " - Pdg Code : " << mcPosV0Dghter->GetTitle() << endl;
3193                 cout << "     - V0 Neg. / Label : " << lblNegV0Dghter
3194                 << " - Pdg Code : " << mcNegV0Dghter->GetTitle() << endl;
3195
3196                 cout << "- Xi daughters - " << endl;
3197                 cout << "     + V0 Pos. mother / Label : " << lblMotherPosV0Dghter
3198                 << " - Pdg Code : " << mcMotherPosV0Dghter->GetTitle() << endl;
3199                 cout << "     - V0 Neg. mother / Label : " << lblMotherNegV0Dghter
3200                 << " - Pdg Code : " << mcMotherNegV0Dghter->GetTitle() << endl;
3201
3202                 cout << "     --  Bach. / Label :" << lblBach
3203                 << " -  Pdg Code : " << mcBach->GetTitle() << endl;
3204
3205                 cout << "- Xi candidate -" << endl;
3206                 cout << "    +  V0 Pos. Gd Mother / Label : " << lblGdMotherPosV0Dghter
3207                 << " - Pdg Code : " << mcGdMotherPosV0Dghter->GetTitle() << endl;
3208                 cout << "    -  V0 Neg. Gd Mother / Label : "  << lblGdMotherNegV0Dghter
3209                 << " - Pdg Code : "<< mcGdMotherNegV0Dghter->GetTitle() << endl;
3210
3211                 cout << "    --  Mother Bach. / Label : " << lblMotherBach
3212                 << " - Pdg Code    : " << mcMotherBach->GetTitle() << endl;
3213                 cout << endl;
3214         }
3215
3216
3217         lmcPt             = mcMotherBach->Pt();
3218         lmcRapCasc        = 0.5*TMath::Log( (mcMotherBach->Energy() + mcMotherBach->Pz()) /
3219                                                      (mcMotherBach->Energy() - mcMotherBach->Pz() +1.e-13) );
3220         lmcEta            = mcMotherBach->Eta();
3221         lmcTransvRadius   = mcBach->R(); // to get the decay point of Xi, = the production vertex of Bachelor ...
3222
3223         TVector3 lmcTVect3Mom( mcMotherBach->Px(), mcMotherBach->Py(), mcMotherBach->Pz() );
3224
3225         lrecoPt           = xiESD->Pt();
3226         lrecoTransvRadius = TMath::Sqrt( xiESD->Xv() * xiESD->Xv() + xiESD->Yv() * xiESD->Yv() );
3227
3228         TVector3 lrecoTVect3Mom( xiESD->Px(), xiESD->Py(), xiESD->Pz() );
3229         lDeltaPhiMcReco   = lmcTVect3Mom.DeltaPhi( lrecoTVect3Mom ) * 180.0/TMath::Pi();
3230
3231         lmcPtPosV0Dghter = mcPosV0Dghter->Pt() ;
3232         lmcPtNegV0Dghter = mcNegV0Dghter->Pt();
3233         lrecoP           = xiESD->P();
3234
3235         Double_t nV0mom[3] = {0. ,0. ,0. };
3236         Double_t pV0mom[3] = {0. ,0. ,0. };
3237
3238         xiESD->GetNPxPyPz(nV0mom[0],nV0mom[1],nV0mom[2]);
3239         xiESD->GetPPxPyPz(pV0mom[0],pV0mom[1],pV0mom[2]);
3240
3241         lV0mom = TMath::Sqrt(TMath::Power(nV0mom[0]+pV0mom[0],2)+TMath::Power(nV0mom[1]+pV0mom[1],2)+TMath::Power(nV0mom[2]+pV0mom[2],2));
3242
3243         Double_t lBachMomX = 0.; Double_t lBachMomY = 0.; Double_t lBachMomZ = 0.;
3244         xiESD->GetBPxPyPz(  lBachMomX,  lBachMomY,  lBachMomZ );
3245         lBachTransvMom   = TMath::Sqrt( lBachMomX*lBachMomX   + lBachMomY*lBachMomY );
3246
3247         lnTrackTransvMom = TMath::Sqrt( nV0mom[0]*nV0mom[0] + nV0mom[1]*nV0mom[1] );
3248         lpTrackTransvMom = TMath::Sqrt( pV0mom[0]*pV0mom[0] + pV0mom[1]*pV0mom[1] );
3249
3250
3251       } else if ( fAnalysisType == "AOD" ) {
3252
3253         const AliAODcascade *xiAOD = lAODevent->GetCascade(iXi);
3254
3255         if (!xiAOD) continue;
3256
3257         // - Step II.1 : Connection to daughter tracks of the current cascade
3258         //-------------
3259
3260         AliAODTrack *pTrackXi    = dynamic_cast<AliAODTrack*>( xiAOD->GetDaughter(0) );
3261         AliAODTrack *nTrackXi    = dynamic_cast<AliAODTrack*>( xiAOD->GetDaughter(1) );
3262         AliAODTrack *bachTrackXi = dynamic_cast<AliAODTrack*>( xiAOD->GetDecayVertexXi()->GetDaughter(0) );
3263         if (!pTrackXi || !nTrackXi || !bachTrackXi ) {
3264           AliWarning("ERROR: Could not retrieve one of the 3 AOD daughter tracks of the cascade ...");
3265           continue;
3266         }
3267
3268         UInt_t lIdxPosXi  = (UInt_t) TMath::Abs( pTrackXi->GetID() );
3269         UInt_t lIdxNegXi  = (UInt_t) TMath::Abs( nTrackXi->GetID() );
3270         UInt_t lBachIdx   = (UInt_t) TMath::Abs( bachTrackXi->GetID() );
3271
3272                 // abs value not needed ; the index should always be positive (!= label ...)
3273
3274
3275         // FIXME : rejection of a double use of a daughter track (nothing but just a crosscheck of what is done in the cascade vertexer)
3276         if(lBachIdx == lIdxNegXi) {
3277                 AliWarning("Pb / Idx(Bach. track) = Idx(Neg. track) ... continue!"); continue;
3278         }
3279         if(lBachIdx == lIdxPosXi) {
3280                 AliWarning("Pb / Idx(Bach. track) = Idx(Pos. track) ... continue!"); continue;
3281         }
3282
3283         lPosTPCClusters   = pTrackXi->GetTPCNcls();
3284         lNegTPCClusters   = nTrackXi->GetTPCNcls();
3285         lBachTPCClusters  = bachTrackXi->GetTPCNcls();
3286
3287
3288         // FIXME : rejection of a poor quality tracks
3289         if (fkQualityCutTPCrefit) {
3290                 // 1 - Poor quality related to TPCrefit
3291           if (!(pTrackXi->IsOn(AliAODTrack::kTPCrefit))) { AliWarning("Pb / V0 Pos. track has no TPCrefit ... continue!"); continue; }
3292           if (!(nTrackXi->IsOn(AliAODTrack::kTPCrefit))) { AliWarning("Pb / V0 Neg. track has no TPCrefit ... continue!"); continue; }
3293           if (!(bachTrackXi->IsOn(AliAODTrack::kTPCrefit))) { AliWarning("Pb / Bach.   track has no TPCrefit ... continue!"); continue; }
3294
3295         }
3296         if (fkQualityCutnTPCcls) {
3297                 // 2 - Poor quality related to TPC clusters
3298                 if(lPosTPCClusters  < fMinnTPCcls) { AliWarning("Pb / V0 Pos. track has less than n TPC clusters ... continue!"); continue; }
3299                 if(lNegTPCClusters  < fMinnTPCcls) { AliWarning("Pb / V0 Neg. track has less than n TPC clusters ... continue!"); continue; }
3300                 if(lBachTPCClusters < fMinnTPCcls) { AliWarning("Pb / Bach.   track has less than n TPC clusters ... continue!"); continue; }
3301         }
3302
3303         etaPos  = pTrackXi->Eta();
3304         etaNeg  = nTrackXi->Eta();
3305         etaBach = bachTrackXi->Eta();
3306
3307         // - Step II.2 : Info over reconstructed cascades
3308         //-------------
3309
3310
3311         if( bachTrackXi->Charge() < 0 ) {
3312                 lInvMassXiMinus = xiAOD->MassXi();
3313
3314                 lInvMassOmegaMinus = xiAOD->MassOmega();
3315         }
3316
3317         if( bachTrackXi->Charge() >  0 ){
3318                 lInvMassXiPlus = xiAOD->MassXi();
3319
3320                 lInvMassOmegaPlus = xiAOD->MassOmega();
3321
3322         }
3323
3324
3325         //lChi2Xi                       = xiAOD->Chi2Xi();
3326         lDcaXiDaughters                 = xiAOD->DcaXiDaughters();
3327         lDcaBachToPrimVertexXi          = xiAOD->DcaBachToPrimVertex();
3328         lXiCosineOfPointingAngle        = xiAOD->CosPointingAngleXi( lBestPrimaryVtxPos[0],
3329                                                                      lBestPrimaryVtxPos[1],
3330                                                                      lBestPrimaryVtxPos[2] );
3331
3332         lPosXi[0] = xiAOD->DecayVertexXiX();
3333         lPosXi[1] = xiAOD->DecayVertexXiY();
3334         lPosXi[2] = xiAOD->DecayVertexXiZ();
3335
3336                                         // This value shouldn't change, whatever the working hyp. is : Xi-, Xi+, Omega-, Omega+
3337         lDcaV0DaughtersXi               = xiAOD->DcaV0Daughters();
3338         // lV0Chi2Xi                    = xiAOD->GetChi2V0();
3339         lV0toXiCosineOfPointingAngle    = xiAOD->CosPointingAngle( lPosXi );
3340
3341         lV0CosineOfPointingAngle        = xiAOD->CosPointingAngle( lBestPrimaryVtxPos );
3342
3343         lPosV0Xi[0] = xiAOD->DecayVertexV0X();
3344         lPosV0Xi[1] = xiAOD->DecayVertexV0Y();
3345         lPosV0Xi[2] = xiAOD->DecayVertexV0Z();
3346
3347
3348         lDcaV0ToPrimVertexXi            = xiAOD->DcaV0ToPrimVertex();
3349
3350         lDcaPosToPrimVertexXi           = xiAOD->DcaPosToPrimVertex();
3351
3352         lDcaNegToPrimVertexXi           = xiAOD->DcaNegToPrimVertex();
3353
3354         lChargeXi = xiAOD->ChargeXi();
3355
3356         // - Step II.3 : PID info
3357
3358         //-------------
3359
3360
3361         // 3.1 - PID Information
3362
3363
3364         // Combined VO-positive-daughter PID
3365
3366         // Combined bachelor PID
3367 /*        AliPID bachPidXi;       bachPidXi.SetPriors(    lPriorsGuessXi    );
3368         AliPID bachPidOmega;    bachPidOmega.SetPriors( lPriorsGuessOmega );
3369
3370         if ( bachTrackXi->IsOn(AliESDtrack::kESDpid) ) {  // Combined PID exists
3371                 Double_t r[10] = {0.}; bachTrackXi->GetESDpid(r);
3372                 bachPidXi.SetProbabilities(r);
3373                 bachPidOmega.SetProbabilities(r);
3374                 // Check if the bachelor track is a pion
3375                     ppionBach = bachPidXi.GetProbability(AliPID::kPion);
3376                 if (ppionBach > bachPidXi.GetProbability(AliPID::kElectron) &&
3377                     ppionBach > bachPidXi.GetProbability(AliPID::kMuon)     &&
3378                     ppionBach > bachPidXi.GetProbability(AliPID::kKaon)     &&
3379                     ppionBach > bachPidXi.GetProbability(AliPID::kProton)   )     lIsBachelorPion = kTRUE;
3380                 // Check if the bachelor track is a kaon
3381                     pkaonBach = bachPidOmega.GetProbability(AliPID::kKaon);
3382                 if (pkaonBach > bachPidOmega.GetProbability(AliPID::kElectron) &&
3383                     pkaonBach > bachPidOmega.GetProbability(AliPID::kMuon)     &&
3384                     pkaonBach > bachPidOmega.GetProbability(AliPID::kPion)     &&
3385                     pkaonBach > bachPidOmega.GetProbability(AliPID::kProton)   )  lIsBachelorKaon = kTRUE;
3386         }// end if bachelor track with existing combined PID
3387 */
3388
3389         // 3.1.B - TPC PID : 4-sigma bands on Bethe-Bloch curve
3390
3391
3392                 // Bachelor
3393         if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( bachTrackXi,AliPID::kKaon)) < 4) lIsBachelorKaonForTPC = kTRUE;
3394         if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( bachTrackXi,AliPID::kPion)) < 4) lIsBachelorPionForTPC = kTRUE;
3395
3396         // Negative V0 daughter
3397         if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( nTrackXi,AliPID::kPion   )) < 4) lIsNegPionForTPC   = kTRUE;
3398         if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( nTrackXi,AliPID::kProton )) < 4) lIsNegProtonForTPC = kTRUE;
3399
3400         // Positive V0 daughter
3401         if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( pTrackXi,AliPID::kPion   )) < 4) lIsPosPionForTPC   = kTRUE;
3402         if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( pTrackXi,AliPID::kProton )) < 4) lIsPosProtonForTPC = kTRUE;
3403
3404 /*
3405         const AliExternalTrackParam *pInnerWallTrackXi    = pTrackXi    ->GetInnerParam(); // Do not use GetTPCInnerWall
3406         const AliExternalTrackParam *nInnerWallTrackXi    = nTrackXi    ->GetInnerParam();
3407         const AliExternalTrackParam *bachInnerWallTrackXi = bachTrackXi ->GetInnerParam();
3408         if(pInnerWallTrackXi && nInnerWallTrackXi && bachInnerWallTrackXi ){
3409
3410                 Double_t pMomInnerWall    = pInnerWallTrackXi   ->GetP();
3411                 Double_t nMomInnerWall    = nInnerWallTrackXi   ->GetP();
3412                 Double_t bachMomInnerWall = bachInnerWallTrackXi->GetP();
3413
3414                 // Bachelor
3415                 if (TMath::Abs(fESDpid->NumberOfSigmasTPC( bachTrackXi,AliPID::kPion)) < 3)                              lIsBachelorPionForTPC = kTRUE;
3416                 if (bachMomInnerWall < 0.350  && TMath::Abs(fESDpid->NumberOfSigmasTPC( bachTrackXi,AliPID::kKaon)) < 5) lIsBachelorKaonForTPC = kTRUE;
3417                 if (bachMomInnerWall > 0.350  && TMath::Abs(fESDpid->NumberOfSigmasTPC( bachTrackXi,AliPID::kKaon)) < 3) lIsBachelorKaonForTPC = kTRUE;
3418
3419                 // Negative V0 daughter
3420                 if (TMath::Abs(fESDpid->NumberOfSigmasTPC( nTrackXi,AliPID::kPion   )) < 3  )                           lIsNegPionForTPC   = kTRUE;
3421                 if (nMomInnerWall < 0.6  && TMath::Abs(fESDpid->NumberOfSigmasTPC( nTrackXi,AliPID::kProton ) ) < 5 )   lIsNegProtonForTPC = kTRUE;
3422                 if (nMomInnerWall > 0.6  && TMath::Abs(fESDpid->NumberOfSigmasTPC( nTrackXi,AliPID::kProton ) ) < 3 )   lIsNegProtonForTPC = kTRUE;
3423
3424                 // Positive V0 daughter
3425                 if (TMath::Abs(fESDpid->NumberOfSigmasTPC( pTrackXi,AliPID::kPion   )) < 3 )                            lIsPosPionForTPC   = kTRUE;
3426                 if (pMomInnerWall < 0.6  && TMath::Abs(fESDpid->NumberOfSigmasTPC( pTrackXi,AliPID::kProton )) < 5)     lIsPosProtonForTPC = kTRUE;
3427                 if (pMomInnerWall > 0.6  && TMath::Abs(fESDpid->NumberOfSigmasTPC( pTrackXi,AliPID::kProton )) < 3)     lIsPosProtonForTPC = kTRUE;
3428
3429         }
3430 */
3431
3432         // 3.2 - PID proba Vs Pt(Bach)
3433         Int_t      lblBachForPID  = (Int_t) TMath::Abs( bachTrackXi->GetLabel() );
3434         AliAODMCParticle* mcBachForPID   = (AliAODMCParticle*) arrayMC->At( lblBachForPID );
3435         lmcPtBach      = mcBachForPID->Pt();
3436
3437         // 3.3 - MC perfect PID
3438
3439         if( mcBachForPID->PdgCode() == -211) lIsBachelorMCPiMinus = kTRUE;
3440         if( mcBachForPID->PdgCode() ==  211) lIsBachelorMCPiPlus  = kTRUE;
3441         if( mcBachForPID->PdgCode() == -321) lIsBachelorMCKMinus  = kTRUE;
3442         if( mcBachForPID->PdgCode() ==  321) lIsBachelorMCKPlus   = kTRUE;
3443
3444
3445
3446         // - Step II.4 : MC association (care : lots of "continue;" below this line)
3447         //-------------
3448
3449         if(fDebug > 5)
3450                 cout    << "MC EventNumber : " << lMCevent->Header()->GetEvent()
3451                         << " / MC event Number in Run : " << lMCevent->Header()->GetEventNrInRun() << endl;
3452
3453         // - Step 4.1 : level of the V0 daughters
3454
3455         Int_t lblPosV0Dghter = (Int_t) TMath::Abs( pTrackXi->GetLabel() );
3456                 // Abs value = needed ! question of quality track association ... (negative label when at least one cluster in the track is from a different particle)
3457         Int_t lblNegV0Dghter = (Int_t) TMath::Abs( nTrackXi->GetLabel() );
3458         AliAODMCParticle* mcPosV0Dghter = (AliAODMCParticle*) arrayMC->At( lblPosV0Dghter );
3459         AliAODMCParticle* mcNegV0Dghter = (AliAODMCParticle*) arrayMC->At( lblNegV0Dghter );
3460
3461
3462         // - Step 4.2 : level of the Xi daughters
3463
3464         Int_t lblMotherPosV0Dghter = mcPosV0Dghter->GetMother();
3465         Int_t lblMotherNegV0Dghter = mcNegV0Dghter->GetMother();
3466
3467                 if( lblMotherPosV0Dghter != lblMotherNegV0Dghter) continue; // same mother
3468                 if( lblMotherPosV0Dghter < 0 ) continue; // this particle is primary, no mother
3469                 if( lblMotherNegV0Dghter < 0 ) continue;
3470
3471
3472                 // mothers = Lambda candidate ... a priori
3473
3474         AliAODMCParticle* mcMotherPosV0Dghter = (AliAODMCParticle*) arrayMC->At( lblMotherPosV0Dghter );
3475         AliAODMCParticle* mcMotherNegV0Dghter = (AliAODMCParticle*) arrayMC->At( lblMotherNegV0Dghter );  // MN: redundant?? already checked that labels are the same...-->same part from stack
3476
3477         Int_t      lblBach  = (Int_t) TMath::Abs( bachTrackXi->GetLabel() );
3478         AliAODMCParticle* mcBach   = (AliAODMCParticle*) arrayMC->At( lblBach );
3479
3480
3481         // - Step 4.3 : level of Xi candidate
3482
3483         Int_t lblGdMotherPosV0Dghter =   mcMotherPosV0Dghter->GetMother() ;
3484         Int_t lblGdMotherNegV0Dghter =   mcMotherNegV0Dghter->GetMother() ;
3485
3486                 if( lblGdMotherPosV0Dghter != lblGdMotherNegV0Dghter ) continue;
3487                 if( lblGdMotherPosV0Dghter < 0 ) continue; // primary lambda ...
3488                 if( lblGdMotherNegV0Dghter < 0 ) continue; // primary lambda ...
3489
3490                 // Gd mothers = Xi candidate ... a priori
3491
3492         AliAODMCParticle* mcGdMotherPosV0Dghter = (AliAODMCParticle*) arrayMC->At( lblGdMotherPosV0Dghter );
3493         AliAODMCParticle* mcGdMotherNegV0Dghter = (AliAODMCParticle*) arrayMC->At( lblGdMotherNegV0Dghter );
3494
3495         Int_t lblMotherBach = (Int_t) TMath::Abs( mcBach->GetMother() );
3496
3497         if( lblMotherBach != lblGdMotherPosV0Dghter ) continue; //same mother for bach and V0 daughters
3498
3499         AliAODMCParticle* mcMotherBach = (AliAODMCParticle*) arrayMC->At( lblMotherBach );
3500
3501                 // Check if cascade is primary
3502
3503         if (!(mcMotherBach->IsPhysicalPrimary())) continue;
3504
3505
3506         kIsNaturalPart = kTRUE;
3507
3508         if ((lblMotherBach>=endOfHijingEvent)&&(endOfHijingEvent!=-1)&&(mcMotherBach->GetMother()<0)) kIsNaturalPart = kFALSE;
3509
3510         // - Step 4.4 : Manage boolean for association
3511
3512         if( mcMotherBach                ->GetPdgCode() ==   3312 &&
3513             mcGdMotherPosV0Dghter       ->GetPdgCode() ==   3312 &&
3514             mcGdMotherNegV0Dghter       ->GetPdgCode() ==   3312)    {  lAssoXiMinus = kTRUE;
3515                                                                         lInvMassLambdaAsCascDghter = xiAOD->MassLambda();
3516                                                                         cascadeMass = 1.321;
3517                                                                       //  Printf("Xi- asso current index = %n ", lblGdMotherPosV0Dghter);
3518                                                                         nAssoXiMinus++; }
3519
3520         else if( mcMotherBach           ->GetPdgCode() ==  -3312 &&
3521             mcGdMotherPosV0Dghter       ->GetPdgCode() ==  -3312 &&
3522             mcGdMotherNegV0Dghter       ->GetPdgCode() ==  -3312)    {  lAssoXiPlus = kTRUE;
3523                                                                         lInvMassLambdaAsCascDghter = xiAOD->MassAntiLambda();
3524                                                                         cascadeMass = 1.321;
3525                                                                         nAssoXiPlus++; }
3526
3527         else if( mcMotherBach           ->GetPdgCode() ==   3334 &&
3528             mcGdMotherPosV0Dghter       ->GetPdgCode() ==   3334 &&
3529             mcGdMotherNegV0Dghter       ->GetPdgCode() ==   3334)    {  lAssoOmegaMinus = kTRUE;
3530                                                                         lInvMassLambdaAsCascDghter = xiAOD->MassLambda();
3531                                                                         cascadeMass = 1.672;
3532                                                                         nAssoOmegaMinus++; }
3533
3534         else if( mcMotherBach           ->GetPdgCode() ==  -3334 &&
3535             mcGdMotherPosV0Dghter       ->GetPdgCode() ==  -3334 &&
3536             mcGdMotherNegV0Dghter       ->GetPdgCode() ==  -3334)    {  lAssoOmegaPlus = kTRUE;
3537                                                                         lInvMassLambdaAsCascDghter = xiAOD->MassAntiLambda();
3538                                                                         cascadeMass = 1.672;
3539                                                                         nAssoOmegaPlus++; }
3540         //-------------
3541
3542         lmcPt             = mcMotherBach->Pt();
3543         lmcRapCasc        = 0.5*TMath::Log( (mcMotherBach->E() + mcMotherBach->Pz()) /
3544                                                      (mcMotherBach->E() - mcMotherBach->Pz() +1.e-13) );
3545         lmcEta            = mcMotherBach->Eta();
3546         Float_t decayCascX = mcBach->Xv();
3547         Float_t decayCascY = mcBach->Yv();
3548         lmcTransvRadius   = TMath::Sqrt(decayCascX*decayCascX+decayCascY*decayCascY); // decay point of Xi, = the production vertex of Bachelor ...
3549
3550         TVector3 lmcTVect3Mom( mcMotherBach->Px(), mcMotherBach->Py(), mcMotherBach->Pz() );
3551
3552         Double_t xiMomX = xiAOD->MomXiX();
3553         Double_t xiMomY = xiAOD->MomXiY();
3554         Double_t xiMomZ = xiAOD->MomXiZ();
3555
3556         lrecoPt           = TMath::Sqrt( xiMomX*xiMomX   + xiMomY*xiMomY );
3557         lrecoTransvRadius = TMath::Sqrt( xiAOD->DecayVertexXiX() * xiAOD->DecayVertexXiX() + xiAOD->DecayVertexXiY() * xiAOD->DecayVertexXiY() );
3558
3559         TVector3 lrecoTVect3Mom( xiMomX, xiMomY, xiMomZ );
3560         lDeltaPhiMcReco   = lmcTVect3Mom.DeltaPhi( lrecoTVect3Mom ) * 180.0/TMath::Pi();
3561
3562         lmcPtPosV0Dghter = mcPosV0Dghter->Pt() ;
3563         lmcPtNegV0Dghter = mcNegV0Dghter->Pt();
3564         lrecoP           = TMath::Sqrt( xiMomX*xiMomX   + xiMomY*xiMomY   + xiMomZ*xiMomZ );;
3565
3566         Double_t lV0momX = xiAOD->MomV0X();
3567         Double_t lV0momY = xiAOD->MomV0Y();
3568         Double_t lV0momZ = xiAOD->MomV0Z();
3569         lV0mom = TMath::Sqrt(TMath::Power(lV0momX,2)+TMath::Power(lV0momY,2)+TMath::Power(lV0momZ,2));
3570
3571
3572         Double_t lBachMomX = xiAOD->MomBachX();
3573         Double_t lBachMomY = xiAOD->MomBachY();
3574         lBachTransvMom  = TMath::Sqrt( lBachMomX*lBachMomX   + lBachMomY*lBachMomY );
3575
3576         Double_t lV0NMomX = xiAOD->MomNegX();
3577         Double_t lV0NMomY = xiAOD->MomNegY();
3578         Double_t lV0PMomX = xiAOD->MomPosX();
3579         Double_t lV0PMomY = xiAOD->MomPosY();
3580
3581         lnTrackTransvMom = TMath::Sqrt( lV0NMomX*lV0NMomX   + lV0NMomY*lV0NMomY );
3582         lpTrackTransvMom = TMath::Sqrt( lV0PMomX*lV0PMomX   + lV0PMomY*lV0PMomY );
3583
3584
3585       }
3586
3587       lXiRadius                       = TMath::Sqrt( lPosXi[0]*lPosXi[0]  +  lPosXi[1]*lPosXi[1] );
3588       lV0RadiusXi                     = TMath::Sqrt( lPosV0Xi[0]*lPosV0Xi[0]  +  lPosV0Xi[1]*lPosV0Xi[1] );
3589
3590
3591       // Cut on pt of the three daughter tracks
3592       if (lBachTransvMom<fMinPtCutOnDaughterTracks) continue;
3593       if (lpTrackTransvMom<fMinPtCutOnDaughterTracks) continue;
3594       if (lnTrackTransvMom<fMinPtCutOnDaughterTracks) continue;
3595
3596       // Cut on pseudorapidity of the three daughter tracks
3597       if (TMath::Abs(etaBach)>fEtaCutOnDaughterTracks) continue;
3598       if (TMath::Abs(etaPos)>fEtaCutOnDaughterTracks) continue;
3599       if (TMath::Abs(etaNeg)>fEtaCutOnDaughterTracks) continue;
3600
3601
3602       // Extra-selection for cascade candidates
3603                 // Towards optimisation of AA selection
3604       if (fkExtraSelections) {
3605
3606         // if(lChi2Xi > 2000) continue; // in AliCascadeVertexer
3607         // if(lV0Chi2Xi > 2000) continue; // in AliV0vertexer
3608
3609         if (lDcaXiDaughters > 0.3) continue; // in AliCascadeVertexer
3610         if (lXiCosineOfPointingAngle < 0.999 ) continue; // in AliCascadeVertexer
3611         if (lDcaV0ToPrimVertexXi < 0.05) continue; // in AliCascadeVertexer
3612         if (lDcaBachToPrimVertexXi < 0.03) continue; // in AliCascadeVertexer
3613 ////      if (TMath::Abs(lInvMassLambdaAsCascDghter-1.11568) > 0.006 ) continue;  // in AliCascadeVertexer
3614
3615         if (lDcaV0DaughtersXi > 1.) continue; // in AliV0vertexer
3616         if (lV0toXiCosineOfPointingAngle < 0.998) continue; // in AliV0vertexer
3617         if (lDcaPosToPrimVertexXi < 0.1) continue; // in AliV0vertexer
3618         if (lDcaNegToPrimVertexXi < 0.1) continue; // in AliV0vertexer
3619
3620
3621           if(lXiRadius < .9) continue; // in AliCascadeVertexer
3622 //        if(lXiRadius > 100) continue; // in AliCascadeVertexer
3623           if(lV0RadiusXi < 0.9) continue; // in AliV0vertexer
3624 //        if(lV0RadiusXi > 100) continue; // in AliV0vertexer
3625
3626       }
3627
3628
3629       // Combined PID TH1s
3630       if( lChargeXi < 0 && lIsBachelorPion )    fHistMassWithCombPIDXiMinus     ->Fill( lInvMassXiMinus    );
3631       if( lChargeXi > 0 && lIsBachelorPion )    fHistMassWithCombPIDXiPlus      ->Fill( lInvMassXiPlus     );
3632       if( lChargeXi < 0 && lIsBachelorKaon )    fHistMassWithCombPIDOmegaMinus  ->Fill( lInvMassOmegaMinus );
3633       if( lChargeXi > 0 && lIsBachelorKaon )    fHistMassWithCombPIDOmegaPlus   ->Fill( lInvMassOmegaPlus  );
3634
3635
3636       if( lChargeXi < 0 )     fHistMassXiMinus        ->Fill( lInvMassXiMinus );
3637       if( lChargeXi > 0 )     fHistMassXiPlus         ->Fill( lInvMassXiPlus );
3638       if( lChargeXi < 0 )     fHistMassOmegaMinus     ->Fill( lInvMassOmegaMinus );
3639       if( lChargeXi > 0 )     fHistMassOmegaPlus      ->Fill( lInvMassOmegaPlus );
3640
3641       if(lIsBachelorPion)   f2dHistPIDprobaPionVsMCPtBach->Fill( lmcPtBach, ppionBach );
3642       if(lIsBachelorKaon)   f2dHistPIDprobaKaonVsMCPtBach->Fill( lmcPtBach, pkaonBach );
3643
3644       if( lChargeXi < 0 && lIsBachelorMCPiMinus )    fHistMassWithMcPIDXiMinus     ->Fill( lInvMassXiMinus );
3645       if( lChargeXi > 0 && lIsBachelorMCPiPlus  )    fHistMassWithMcPIDXiPlus      ->Fill( lInvMassXiPlus );
3646       if( lChargeXi < 0 && lIsBachelorMCKMinus  )    fHistMassWithMcPIDOmegaMinus  ->Fill( lInvMassOmegaMinus );
3647       if( lChargeXi > 0 && lIsBachelorMCKPlus   )    fHistMassWithMcPIDOmegaPlus   ->Fill( lInvMassOmegaPlus );
3648
3649
3650
3651       if(!lAssoXiMinus && !lAssoXiPlus && !lAssoOmegaMinus && !lAssoOmegaPlus) continue; // no association, skip the rest of the code
3652
3653       // - Fill pt histos. xim only!!
3654       if (lAssoXiMinus) {
3655         fHistPtRecBachXiMinus->Fill(lBachTransvMom);
3656         fHistPtRecMesDghterXiMinus->Fill(lpTrackTransvMom);
3657         fHistPtRecBarDghterXiMinus->Fill(lnTrackTransvMom);
3658       }
3659
3660         // Calculate proper time for cascade (reconstructed)
3661
3662         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));
3663         if (lrecoP!=0)         lctau = lctau*cascadeMass/lrecoP;
3664         else lctau = -1.;
3665
3666
3667         // Calculate proper time for Lambda (reconstructed)
3668         Float_t lambdaMass = 1.115683; // PDG mass
3669         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));
3670         Float_t lctauV0 = -1.;
3671         if (lV0mom!=0) lctauV0 = distV0Xi*lambdaMass/lV0mom;
3672
3673         Float_t distTV0Xi =  TMath::Sqrt(TMath::Power((lPosV0Xi[0]-lPosXi[0]),2)+TMath::Power((lPosV0Xi[1]-lPosXi[1]),2));
3674
3675
3676         // - Histos for the cascade candidates associated with MC
3677
3678         if( lChargeXi < 0 && lAssoXiMinus){
3679                 fHistAsMCMassXiMinus          ->Fill( lInvMassXiMinus  );
3680                 if(lIsBachelorPion)     f2dHistAsMCandCombPIDGenPtVsGenYXiMinus->Fill( lmcPt, lmcRapCasc );
3681                 f2dHistAsMCGenPtVsGenYXiMinus ->Fill( lmcPt, lmcRapCasc);
3682                 fHistAsMCGenEtaXiMinus        ->Fill( lmcEta           );
3683                 f2dHistAsMCResPtXiMinus       ->Fill( lmcPt,           (lrecoPt - lmcPt)/ lmcPt );
3684                 f2dHistAsMCResRXiMinus        ->Fill( lmcTransvRadius, (lrecoTransvRadius - lmcTransvRadius)/ lmcTransvRadius    );
3685                 f2dHistAsMCResPhiXiMinus      ->Fill( lmcPt, lDeltaPhiMcReco );
3686                 f2dHistAsMCptProtonMCptXiMinus->Fill(lmcPt,lmcPtPosV0Dghter);
3687                 fHistV0CosineOfPointingAnglevsPtXi->Fill(lmcPt,lV0CosineOfPointingAngle);
3688         }
3689
3690         else if( lChargeXi > 0 && lAssoXiPlus){
3691                 fHistAsMCMassXiPlus           ->Fill( lInvMassXiPlus   );
3692                 if(lIsBachelorPion)     f2dHistAsMCandCombPIDGenPtVsGenYXiPlus->Fill( lmcPt, lmcRapCasc );
3693                 f2dHistAsMCGenPtVsGenYXiPlus  ->Fill( lmcPt, lmcRapCasc);
3694                 fHistAsMCGenEtaXiPlus         ->Fill( lmcEta           );
3695                 f2dHistAsMCResPtXiPlus        ->Fill( lmcPt,           (lrecoPt - lmcPt)/ lmcPt );
3696                 f2dHistAsMCResRXiPlus         ->Fill( lmcTransvRadius, (lrecoTransvRadius - lmcTransvRadius)/ lmcTransvRadius    );
3697                 f2dHistAsMCResPhiXiPlus       ->Fill( lmcPt, lDeltaPhiMcReco );
3698                 f2dHistAsMCptAntiprotonMCptXiPlus->Fill(lmcPt,lmcPtNegV0Dghter);
3699                 fHistV0CosineOfPointingAnglevsPtXi->Fill(lmcPt,lV0CosineOfPointingAngle);
3700         }
3701
3702         else if( lChargeXi < 0 && lAssoOmegaMinus){
3703                 fHistAsMCMassOmegaMinus          ->Fill( lInvMassOmegaMinus );
3704                 if(lIsBachelorKaon)     f2dHistAsMCandCombPIDGenPtVsGenYOmegaMinus->Fill( lmcPt, lmcRapCasc );
3705                 f2dHistAsMCGenPtVsGenYOmegaMinus ->Fill( lmcPt, lmcRapCasc  );
3706                 fHistAsMCGenEtaOmegaMinus        ->Fill( lmcEta             );
3707                 f2dHistAsMCResPtOmegaMinus       ->Fill( lmcPt,           (lrecoPt - lmcPt)/ lmcPt );
3708                 f2dHistAsMCResROmegaMinus        ->Fill( lmcTransvRadius, (lrecoTransvRadius - lmcTransvRadius)/ lmcTransvRadius    );
3709                 f2dHistAsMCResPhiOmegaMinus      ->Fill( lmcPt, lDeltaPhiMcReco );
3710                 f2dHistAsMCptProtonMCptOmegaMinus->Fill(lmcPt,lmcPtPosV0Dghter);
3711                 fHistV0CosineOfPointingAnglevsPtOmega->Fill(lmcPt,lV0CosineOfPointingAngle);
3712         }
3713
3714         else if( lChargeXi > 0 && lAssoOmegaPlus){
3715                 fHistAsMCMassOmegaPlus           ->Fill( lInvMassOmegaPlus );
3716                 if(lIsBachelorKaon)     f2dHistAsMCandCombPIDGenPtVsGenYOmegaPlus->Fill( lmcPt, lmcRapCasc );
3717                 f2dHistAsMCGenPtVsGenYOmegaPlus  ->Fill( lmcPt, lmcRapCasc   );
3718                 fHistAsMCGenEtaOmegaPlus         ->Fill( lmcEta            );
3719                 f2dHistAsMCResPtOmegaPlus        ->Fill( lmcPt,           (lrecoPt - lmcPt)/ lmcPt );
3720                 f2dHistAsMCResROmegaPlus         ->Fill( lmcTransvRadius, (lrecoTransvRadius - lmcTransvRadius)/ lmcTransvRadius    );
3721                 f2dHistAsMCResPhiOmegaPlus       ->Fill( lmcPt, lDeltaPhiMcReco );
3722                 f2dHistAsMCptAntiprotonMCptOmegaPlus->Fill(lmcPt,lmcPtNegV0Dghter);
3723                 fHistV0CosineOfPointingAnglevsPtOmega->Fill(lmcPt,lV0CosineOfPointingAngle);
3724         }
3725
3726         fHistV0toXiCosineOfPointingAngle->Fill(lV0toXiCosineOfPointingAngle);
3727
3728         // - Step 6 : Containers = Cascade cuts + PID
3729         //-------------
3730         // 6.3 - Filling the AliCFContainer (optimisation of topological selections + systematics)
3731         Double_t lContainerCutVars[22] = {0.0};
3732
3733         lContainerCutVars[0]  = lDcaXiDaughters;
3734         lContainerCutVars[1]  = lDcaBachToPrimVertexXi;
3735         lContainerCutVars[2]  = lXiCosineOfPointingAngle;
3736         lContainerCutVars[3]  = lXiRadius;
3737         lContainerCutVars[4]  = lInvMassLambdaAsCascDghter;
3738         lContainerCutVars[5]  = lDcaV0DaughtersXi;
3739         lContainerCutVars[6]  = lV0CosineOfPointingAngle;
3740         lContainerCutVars[7]  = lV0RadiusXi;
3741         lContainerCutVars[8]  = lDcaV0ToPrimVertexXi;
3742         lContainerCutVars[9]  = lDcaPosToPrimVertexXi;
3743         lContainerCutVars[10] = lDcaNegToPrimVertexXi;
3744
3745         lContainerCutVars[13] = lmcPt;
3746
3747         if (kIsNaturalPart) lContainerCutVars[16] = 0.;
3748         else lContainerCutVars[16] = 1.;
3749         lContainerCutVars[17] = lcentrality;
3750         lContainerCutVars[18] = lPrimaryTrackMultiplicity;
3751         lContainerCutVars[19] = lctau;
3752         lContainerCutVars[20] = lctauV0;
3753         lContainerCutVars[21] = distTV0Xi;
3754
3755         // All cases should be covered below
3756         if( lChargeXi < 0 && lAssoXiMinus    ) {
3757                 lContainerCutVars[11] = lInvMassXiMinus;
3758                 lContainerCutVars[12] = lInvMassOmegaMinus;//1.63;
3759                 lContainerCutVars[14] = lmcRapCasc;
3760                 lContainerCutVars[15] = -1.;
3761                 if ( lIsBachelorPionForTPC   && lIsPosProtonForTPC    && lIsNegPionForTPC ) {
3762                   fCFContAsCascadeCuts->Fill(lContainerCutVars,0); // for Xi-
3763                   fHistEtaBachXiM->Fill(etaBach);
3764                   fHistEtaPosXiM->Fill(etaPos);
3765                   fHistEtaNegXiM->Fill(etaNeg);
3766                 }
3767         }
3768         if( lChargeXi > 0 && lAssoXiPlus    ) {
3769                 lContainerCutVars[11] = lInvMassXiPlus;
3770                 lContainerCutVars[12] = lInvMassOmegaPlus;//1.26;
3771                 lContainerCutVars[14] = lmcRapCasc;
3772                 lContainerCutVars[15] = -1.;
3773                 if ( lIsBachelorPionForTPC   && lIsNegProtonForTPC    && lIsPosPionForTPC )
3774                   fCFContAsCascadeCuts->Fill(lContainerCutVars,1); // for Xi+
3775         }
3776         if( lChargeXi < 0 && lAssoOmegaMinus ) {
3777                 lContainerCutVars[11] = lInvMassXiMinus;//1.63;
3778                 lContainerCutVars[12] = lInvMassOmegaMinus;
3779                 lContainerCutVars[14] = -1.;
3780                 lContainerCutVars[15] = lmcRapCasc;
3781                 if ( lIsBachelorKaonForTPC   && lIsPosProtonForTPC    && lIsNegPionForTPC )
3782                   fCFContAsCascadeCuts->Fill(lContainerCutVars,2); // for Omega-
3783         }
3784         if( lChargeXi > 0 && lAssoOmegaPlus  ) {
3785                 lContainerCutVars[11] = lInvMassXiPlus;//1.26;
3786                 lContainerCutVars[12] = lInvMassOmegaPlus;
3787                 lContainerCutVars[14] = -1.;
3788                 lContainerCutVars[15] = lmcRapCasc;
3789                 if ( lIsBachelorKaonForTPC   && lIsNegProtonForTPC    && lIsPosPionForTPC )
3790                   fCFContAsCascadeCuts->Fill(lContainerCutVars,3); // for Omega+
3791         }
3792
3793
3794         // 6.4 - Filling the AliCFContainers related to PID
3795
3796         Double_t lContainerPIDVars[4] = {0.0};
3797
3798
3799         // Xi Minus
3800         if( lChargeXi < 0 && lAssoXiMinus ) {
3801                 lContainerPIDVars[0] = lmcPt              ;
3802                 lContainerPIDVars[1] = lInvMassXiMinus    ;
3803                 lContainerPIDVars[2] = lmcRapCasc         ;
3804                 lContainerPIDVars[3] = lcentrality;
3805
3806                 // No PID
3807                         fCFContCascadePIDAsXiMinus->Fill(lContainerPIDVars, 0); // No PID
3808                 // TPC PID
3809                 if( lIsBachelorPionForTPC  )
3810                         fCFContCascadePIDAsXiMinus->Fill(lContainerPIDVars, 1); // TPC PID / 4-#sigma cut on Bachelor track
3811
3812                 if( lIsBachelorPionForTPC &&
3813                     lIsPosProtonForTPC     )
3814                         fCFContCascadePIDAsXiMinus->Fill(lContainerPIDVars, 2); // TPC PID / 4-#sigma cut on Bachelor+Baryon tracks
3815
3816                 if( lIsBachelorPionForTPC &&
3817                     lIsPosProtonForTPC    &&
3818                     lIsNegPionForTPC       )
3819                         fCFContCascadePIDAsXiMinus->Fill(lContainerPIDVars, 3); // TPC PID / 4-#sigma cut on Bachelor+Baryon+Meson tracks
3820
3821                 // Combined PID
3822                 if( lIsBachelorPion        )
3823                         fCFContCascadePIDAsXiMinus->Fill(lContainerPIDVars, 4); // Comb. PID / Bachelor
3824
3825                 if( lIsBachelorPion       &&
3826                     lIsPosInXiProton    )
3827                         fCFContCascadePIDAsXiMinus->Fill(lContainerPIDVars, 5); // Comb. PID / Bachelor+Baryon
3828
3829                 if( lIsBachelorPion     &&
3830                     lIsPosInXiProton &&
3831                     lIsNegInXiPion    )
3832                         fCFContCascadePIDAsXiMinus->Fill(lContainerPIDVars, 6); // Comb. PID / Bachelor+Baryon+Meson
3833         }
3834
3835         lContainerPIDVars[0] = 0.; lContainerPIDVars[1] = 0.; lContainerPIDVars[2] = 0.; lContainerPIDVars[3] = 0.;
3836
3837         // Xi Plus
3838         if( lChargeXi > 0 && lAssoXiPlus ) {
3839                 lContainerPIDVars[0] = lmcPt         ;
3840                 lContainerPIDVars[1] = lInvMassXiPlus;
3841                 lContainerPIDVars[2] = lmcRapCasc    ;
3842                 lContainerPIDVars[3] = lcentrality   ;
3843
3844                 // No PID
3845                         fCFContCascadePIDAsXiPlus->Fill(lContainerPIDVars, 0); // No PID
3846                 // TPC PID
3847                 if( lIsBachelorPionForTPC  )
3848                         fCFContCascadePIDAsXiPlus->Fill(lContainerPIDVars, 1); // TPC PID / 4-#sigma cut on Bachelor track
3849
3850                 if( lIsBachelorPionForTPC &&
3851                     lIsNegProtonForTPC     )
3852                         fCFContCascadePIDAsXiPlus->Fill(lContainerPIDVars, 2); // TPC PID / 4-#sigma cut on Bachelor+Baryon tracks
3853
3854                 if( lIsBachelorPionForTPC &&
3855                     lIsNegProtonForTPC    &&
3856                     lIsPosPionForTPC       )
3857                         fCFContCascadePIDAsXiPlus->Fill(lContainerPIDVars, 3); // TPC PID / 4-#sigma cut on Bachelor+Baryon+Meson tracks
3858
3859                 // Combined PID
3860                 if( lIsBachelorPion        )
3861                         fCFContCascadePIDAsXiPlus->Fill(lContainerPIDVars, 4); // Comb. PID / Bachelor
3862
3863                 if( lIsBachelorPion       &&
3864                     lIsNegInXiProton    )
3865                         fCFContCascadePIDAsXiPlus->Fill(lContainerPIDVars, 5); // Comb. PID / Bachelor+Baryon
3866
3867                 if( lIsBachelorPion     &&
3868                     lIsNegInXiProton &&
3869                     lIsPosInXiPion    )
3870                         fCFContCascadePIDAsXiPlus->Fill(lContainerPIDVars, 6); // Comb. PID / Bachelor+Baryon+Meson
3871         }
3872
3873         lContainerPIDVars[0] = 0.; lContainerPIDVars[1] = 0.; lContainerPIDVars[2] = 0.; lContainerPIDVars[3] = 0.;
3874
3875         // Omega Minus
3876         if( lChargeXi < 0 && lAssoOmegaMinus ) {
3877                 lContainerPIDVars[0] = lmcPt              ;
3878                 lContainerPIDVars[1] = lInvMassOmegaMinus ;
3879                 lContainerPIDVars[2] = lmcRapCasc         ;
3880                 lContainerPIDVars[3] = lcentrality;
3881
3882                 // No PID
3883                         fCFContCascadePIDAsOmegaMinus->Fill(lContainerPIDVars, 0); // No PID
3884                 // TPC PID
3885                 if( lIsBachelorKaonForTPC  )
3886                         fCFContCascadePIDAsOmegaMinus->Fill(lContainerPIDVars, 1); // TPC PID / 4-#sigma cut on Bachelor track
3887
3888                 if( lIsBachelorKaonForTPC &&
3889                     lIsPosProtonForTPC     )
3890                         fCFContCascadePIDAsOmegaMinus->Fill(lContainerPIDVars, 2); // TPC PID / 4-#sigma cut on Bachelor+Baryon tracks
3891
3892                 if( lIsBachelorKaonForTPC &&
3893                     lIsPosProtonForTPC    &&
3894                     lIsNegPionForTPC       )
3895                         fCFContCascadePIDAsOmegaMinus->Fill(lContainerPIDVars, 3); // TPC PID / 4-#sigma cut on Bachelor+Baryon+Meson tracks
3896
3897                 // Combined PID
3898                 if( lIsBachelorKaon        )
3899                         fCFContCascadePIDAsOmegaMinus->Fill(lContainerPIDVars, 4); // Comb. PID / Bachelor
3900
3901                 if( lIsBachelorKaon       &&
3902                     lIsPosInOmegaProton    )
3903                         fCFContCascadePIDAsOmegaMinus->Fill(lContainerPIDVars, 5); // Comb. PID / Bachelor+Baryon
3904
3905                 if(lIsBachelorKaon     &&
3906                    lIsPosInOmegaProton &&
3907                    lIsNegInOmegaPion    )
3908                         fCFContCascadePIDAsOmegaMinus->Fill(lContainerPIDVars, 6); // Comb. PID / Bachelor+Baryon+Meson
3909         }
3910
3911         lContainerPIDVars[0] = 0.; lContainerPIDVars[1] = 0.; lContainerPIDVars[2] = 0.; lContainerPIDVars[3] = 0.;
3912
3913         // Omega Plus
3914         if( lChargeXi > 0 && lAssoOmegaPlus) {
3915                 lContainerPIDVars[0] = lmcPt              ;
3916                 lContainerPIDVars[1] = lInvMassOmegaPlus  ;
3917                 lContainerPIDVars[2] = lmcRapCasc         ;
3918                 lContainerPIDVars[3] = lcentrality;
3919
3920                 // No PID
3921                         fCFContCascadePIDAsOmegaPlus->Fill(lContainerPIDVars, 0); // No PID
3922                 // TPC PID
3923                 if( lIsBachelorKaonForTPC  )
3924                         fCFContCascadePIDAsOmegaPlus->Fill(lContainerPIDVars, 1); // TPC PID / 4-#sigma cut on Bachelor track
3925
3926                 if( lIsBachelorKaonForTPC &&
3927                     lIsNegProtonForTPC     )
3928                         fCFContCascadePIDAsOmegaPlus->Fill(lContainerPIDVars, 2); // TPC PID / 4-#sigma cut on Bachelor+Baryon tracks
3929
3930                 if( lIsBachelorKaonForTPC &&
3931                     lIsNegProtonForTPC    &&
3932                     lIsPosPionForTPC       )
3933                         fCFContCascadePIDAsOmegaPlus->Fill(lContainerPIDVars, 3); // TPC PID / 4-#sigma cut on Bachelor+Baryon+Meson tracks
3934
3935                 // Combined PID
3936                 if( lIsBachelorKaon        )
3937                         fCFContCascadePIDAsOmegaPlus->Fill(lContainerPIDVars, 4); // Comb. PID / Bachelor
3938
3939                 if( lIsBachelorKaon       &&
3940                     lIsNegInOmegaProton    )
3941                         fCFContCascadePIDAsOmegaPlus->Fill(lContainerPIDVars, 5); // Comb. PID / Bachelor+Baryon
3942
3943                 if( lIsBachelorKaon     &&
3944                     lIsNegInOmegaProton &&
3945                     lIsPosInOmegaPion    )
3946                         fCFContCascadePIDAsOmegaPlus->Fill(lContainerPIDVars, 6); // Comb. PID / Bachelor+Baryon+Meson
3947
3948       }
3949
3950 }// End of loop over reconstructed cascades
3951
3952 fHistnAssoXiMinus->Fill(nAssoXiMinus);
3953 fHistnAssoXiPlus->Fill(nAssoXiPlus);
3954 fHistnAssoOmegaMinus->Fill(nAssoOmegaMinus);
3955 fHistnAssoOmegaPlus->Fill(nAssoOmegaPlus);
3956 //Printf("N asso Xi- = %n ", nAssoXiMinus);
3957
3958   // Post output data.
3959  PostData(1, fListHistCascade);
3960  PostData(2, fCFContCascadePIDAsXiMinus);
3961  PostData(3, fCFContCascadePIDAsXiPlus);
3962  PostData(4, fCFContCascadePIDAsOmegaMinus);
3963  PostData(5, fCFContCascadePIDAsOmegaPlus);
3964  PostData(6, fCFContAsCascadeCuts);
3965
3966 }
3967
3968
3969 //________________________________________________________________________
3970 void AliAnalysisTaskCheckPerformanceCascadePbPb::Terminate(Option_t *) {
3971   // Draw result to the screen
3972   // Called once at the end of the query
3973
3974   TList *cRetrievedList = 0x0;
3975   cRetrievedList = (TList*)GetOutputData(1);
3976   if(!cRetrievedList) {
3977         Printf("ERROR - AliAnalysisTaskCheckPerformanceCascadePbPb : ouput data container list not available\n");
3978         return;
3979   }
3980
3981   fHistMCTrackMultiplicity = dynamic_cast<TH1F*> (  cRetrievedList->FindObject("fHistMCTrackMultiplicity")  );
3982   if (!fHistMCTrackMultiplicity) {
3983     Printf("ERROR - AliAnalysisTaskCheckPerformanceCascadePbPb : fHistMCTrackMultiplicity not available");
3984     return;
3985   }
3986
3987
3988   TCanvas *canCheckPerformanceCascade = new TCanvas("CheckPerformanceCascadePbPb","Multiplicity",10,10,510,510);
3989   canCheckPerformanceCascade->cd(1)->SetLogy();
3990
3991   fHistMCTrackMultiplicity->SetMarkerStyle(22);
3992   fHistMCTrackMultiplicity->DrawCopy("E");
3993
3994 }