PWG2/SPECTRA/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.
  18 //            Use with AOD tree = under development
  19 //            Origin   : AliAnalysisTaskCheckPerformanceCascade class by A. Maire Nov2010, antonin.maire@ires.in2p3.fr
  20 //            Modified for PbPb analysis: M. Nicassio Feb2011, maria.nicassio@ba.infn.it:
  21 //                        - physics selection moved to the runProof.C macro
  22 //                        - added centrality selection
  23 //                        - added new histograms
  24 //                        - modified binning of some histograms and containers
  25 //                        - flag to enable CF container usage
  26 //                        - protection in the destructor for CAF usage
  27 //                        - flag for acceptance cut in the MC part
  28 //                        - in the MC particle selection IsPhysicalPrimary added and number of particles taken as appropriate for HIJING
  29 //                          (however for cascades one gets the same if runs on Nprimaries in the stack and does not require IsPhysicalPrimary)
  30 //                        - number of tracklets from AOD also
  31 //                        - automatic settings for PID (July 2011)
  32 //                        - added possibility to select either injected cascades or HIJING cascades (May 2011)
  33 //-----------------------------------------------------------------
  34
  35
  36 #include <Riostream.h>
  37
  38 #include "TList.h"
  39 #include "TFile.h"
  40 #include "TH1F.h"
  41 #include "TH2F.h"
  42 #include "TH3F.h"
  43 #include "TVector3.h"
  44 #include "TCanvas.h"
  45 #include "TParticle.h"
  46 #include "TMath.h"
  47
  48 #include "AliLog.h"
  49 #include "AliHeader.h"
  50 #include "AliMCEvent.h"
  51 #include "AliStack.h"
  52 #include "AliMultiplicity.h"
  53 #include "AliInputEventHandler.h"
  54 #include "AliAnalysisManager.h"
  55
  56 #include "AliCentrality.h"
  57
  58 #include "AliCFContainer.h"
  59
  60 #include "AliESDVZERO.h"
  61
  62 #include "AliGenEventHeader.h"
  63 #include "AliGenCocktailEventHeader.h"
  64 #include "AliGenHijingEventHeader.h"
  65 #include "AliESDtrackCuts.h"
  66 #include "AliPIDResponse.h"
  67 //   #include "AliV0vertexer.h"
  68 //   #include "AliCascadeVertexer.h"
  69 #include "AliESDEvent.h"
  70 #include "AliESDcascade.h"
  71 #include "AliAODEvent.h"
  72 #include "AliAnalysisTaskCheckPerformanceCascadePbPb.h"
  73
  74 ClassImp(AliAnalysisTaskCheckPerformanceCascadePbPb)
  75
  76
  77
  78      //_____Dummy constructor________________________________________________________________
  79 AliAnalysisTaskCheckPerformanceCascadePbPb::AliAnalysisTaskCheckPerformanceCascadePbPb()
  80 : AliAnalysisTaskSE(), // <- take care to AliAnalysisTask( empty )
  81   fDebugCascade(0), fAnalysisType("ESD"), fCollidingSystems(0), fESDtrackCuts(0), /*fPaveTextBookKeeping(0),*/
  82
  83     fPIDResponse                   (0),
  84     fkRerunV0CascVertexers         (0),
  85     fkQualityCutZprimVtxPos        (kTRUE),
  86     fkRejectEventPileUp            (kTRUE),
  87     fkQualityCutNoTPConlyPrimVtx   (kTRUE),
  88     fkQualityCutTPCrefit           (kTRUE),
  89     fkQualityCut80TPCcls           (kTRUE),
  90     fkIsDataRecoWith1PadTPCCluster (kTRUE),
  91     fkExtraSelections              (0),
  92     fUseCFCont                     (0),
  93     fCentrLowLim(0),    fCentrUpLim(0), fCentrEstimator(0),
  94     fVtxRange                      (0),
  95     fApplyAccCut                   (0),
  96
  97
  98         // - Cascade part initialisation
  99     fListHistCascade(0),
 100
 101     // Events in centrality bins
 102     fHistEvtsInCentralityBinsvsNtracks(0),
 103
 104     // Cascade multiplicity histos
 105
 106     fHistnXiPlusPerEvTot(0),
 107     fHistnXiMinusPerEvTot(0),
 108     fHistnOmegaPlusPerEvTot(0),
 109     fHistnOmegaMinusPerEvTot(0),
 110
 111     fHistnXiPlusPerEv(0),
 112     fHistnXiMinusPerEv(0),
 113     fHistnOmegaPlusPerEv(0),
 114     fHistnOmegaMinusPerEv(0),
 115
 116     fHistnAssoXiMinus(0),
 117     fHistnAssoXiPlus(0),
 118     fHistnAssoOmegaMinus(0),
 119     fHistnAssoOmegaPlus(0),
 120
 121
 122     fHistMCTrackMultiplicity(0),
 123        // - Resolution of the multiplicity estimator
 124     f2dHistRecoMultVsMCMult(0),
 125
 126     fHistEtaGenProton(0),
 127     fHistEtaGenAntiProton(0),
 128
 129    // Xi-
 130    fHistEtaGenCascXiMinus(0),
 131    f3dHistGenPtVsGenYGenvsCentXiMinus(0),
 132    f3dHistGenPtVsGenYGenvsNtracksXiMinus(0),
 133
 134
 135     fHistThetaGenCascXiMinus(0),
 136     f2dHistGenPtVsGenYFdblXiMinus(0),
 137
 138     fHistThetaLambdaXiMinus(0),
 139     fHistThetaBachXiMinus(0),
 140
 141     fHistThetaMesDghterXiMinus(0),
 142     fHistThetaBarDghterXiMinus(0),
 143
 144     fHistPtBachXiMinus(0),
 145     fHistPtMesDghterXiMinus(0),
 146     fHistPtBarDghterXiMinus(0),
 147
 148
 149    // Xi+
 150    fHistEtaGenCascXiPlus(0),
 151    f3dHistGenPtVsGenYGenvsCentXiPlus(0),
 152    f3dHistGenPtVsGenYGenvsNtracksXiPlus(0),
 153
 154
 155     fHistThetaGenCascXiPlus(0),
 156     f2dHistGenPtVsGenYFdblXiPlus(0),
 157
 158     fHistThetaLambdaXiPlus(0),
 159     fHistThetaBachXiPlus(0),
 160
 161     fHistThetaMesDghterXiPlus(0),
 162     fHistThetaBarDghterXiPlus(0),
 163
 164     fHistPtBachXiPlus(0),
 165     fHistPtMesDghterXiPlus(0),
 166     fHistPtBarDghterXiPlus(0),
 167
 168    // Omega-
 169    fHistEtaGenCascOmegaMinus(0),
 170    f3dHistGenPtVsGenYGenvsCentOmegaMinus(0),
 171    f3dHistGenPtVsGenYGenvsNtracksOmegaMinus(0),
 172
 173
 174     fHistThetaGenCascOmegaMinus(0),
 175     f2dHistGenPtVsGenYFdblOmegaMinus(0),
 176
 177     fHistThetaLambdaOmegaMinus(0),
 178     fHistThetaBachOmegaMinus(0),
 179
 180     fHistThetaMesDghterOmegaMinus(0),
 181     fHistThetaBarDghterOmegaMinus(0),
 182
 183     fHistPtBachOmegaMinus(0),
 184     fHistPtMesDghterOmegaMinus(0),
 185     fHistPtBarDghterOmegaMinus(0),
 186
 187
 188    // Omega+
 189    fHistEtaGenCascOmegaPlus(0),
 190    f3dHistGenPtVsGenYGenvsCentOmegaPlus(0),
 191    f3dHistGenPtVsGenYGenvsNtracksOmegaPlus(0),
 192
 193
 194     fHistThetaGenCascOmegaPlus(0),
 195     f2dHistGenPtVsGenYFdblOmegaPlus(0),
 196
 197     fHistThetaLambdaOmegaPlus(0),
 198     fHistThetaBachOmegaPlus(0),
 199
 200     fHistThetaMesDghterOmegaPlus(0),
 201     fHistThetaBarDghterOmegaPlus(0),
 202
 203     fHistPtBachOmegaPlus(0),
 204     fHistPtMesDghterOmegaPlus(0),
 205     fHistPtBarDghterOmegaPlus(0),
 206
 207 // Part 2 - Association to MC
 208
 209     fHistMassXiMinus(0),
 210     fHistMassXiPlus(0),
 211     fHistMassOmegaMinus(0),
 212     fHistMassOmegaPlus(0),
 213
 214         // - Effective mass histos with combined PID
 215     fHistMassWithCombPIDXiMinus(0), fHistMassWithCombPIDXiPlus(0),
 216     fHistMassWithCombPIDOmegaMinus(0), fHistMassWithCombPIDOmegaPlus(0),
 217
 218         // - PID Probability versus MC Pt(bachelor track)
 219     f2dHistPIDprobaKaonVsMCPtBach(0), f2dHistPIDprobaPionVsMCPtBach(0),
 220
 221         // - Effective mass histos with perfect MC PID on the bachelor
 222     fHistMassWithMcPIDXiMinus(0), fHistMassWithMcPIDXiPlus(0),
 223     fHistMassWithMcPIDOmegaMinus(0), fHistMassWithMcPIDOmegaPlus(0),
 224
 225
 226         // - Effective mass histos for the cascade candidates associated with MC
 227     fHistAsMCMassXiMinus(0),
 228     fHistAsMCMassXiPlus(0),
 229     fHistAsMCMassOmegaMinus(0),
 230     fHistAsMCMassOmegaPlus(0),
 231
 232         // - Generated Pt Vs generated y, for the cascade candidates associated with MC + Info Comb. PID
 233     f2dHistAsMCandCombPIDGenPtVsGenYXiMinus(0),
 234     f2dHistAsMCandCombPIDGenPtVsGenYXiPlus(0),
 235     f2dHistAsMCandCombPIDGenPtVsGenYOmegaMinus(0),
 236     f2dHistAsMCandCombPIDGenPtVsGenYOmegaPlus(0),
 237
 238         // - Generated Pt Vs generated y, for the cascade candidates associated with MC
 239     f2dHistAsMCGenPtVsGenYXiMinus(0),
 240     f2dHistAsMCGenPtVsGenYXiPlus(0),
 241     f2dHistAsMCGenPtVsGenYOmegaMinus(0),
 242     f2dHistAsMCGenPtVsGenYOmegaPlus(0),
 243
 244         // - Generated Eta of the the cascade candidates associated with MC
 245     fHistAsMCGenEtaXiMinus(0),
 246     fHistAsMCGenEtaXiPlus(0),
 247     fHistAsMCGenEtaOmegaMinus(0),
 248     fHistAsMCGenEtaOmegaPlus(0),
 249
 250         // - Resolution in Pt as function of generated Pt
 251     f2dHistAsMCResPtXiMinus(0),
 252     f2dHistAsMCResPtXiPlus(0),
 253     f2dHistAsMCResPtOmegaMinus(0),
 254     f2dHistAsMCResPtOmegaPlus(0),
 255
 256         // - Resolution in R(2D) as function of generated R
 257     f2dHistAsMCResRXiMinus(0),
 258     f2dHistAsMCResRXiPlus(0),
 259     f2dHistAsMCResROmegaMinus(0),
 260     f2dHistAsMCResROmegaPlus(0),
 261
 262         // - Resolution in phi as function of generated Pt
 263     f2dHistAsMCResPhiXiMinus(0),
 264     f2dHistAsMCResPhiXiPlus(0),
 265     f2dHistAsMCResPhiOmegaMinus(0),
 266     f2dHistAsMCResPhiOmegaPlus(0),
 267
 268     fCFContCascadePIDAsXiMinus(0),
 269     fCFContCascadePIDAsXiPlus(0),
 270     fCFContCascadePIDAsOmegaMinus(0),
 271     fCFContCascadePIDAsOmegaPlus(0),
 272
 273     fCFContAsCascadeCuts(0),
 274
 275     fV0Ampl             (0)
 276
 277 {
 278 // Dummy constructor
 279         for(Int_t iAlephIdx   = 0; iAlephIdx   < 5; iAlephIdx++   ) { fAlephParameters [iAlephIdx]    = -1.; }
 280         for(Int_t iV0selIdx   = 0; iV0selIdx   < 7; iV0selIdx++   ) { fV0Sels          [iV0selIdx   ] = -1.; }
 281         for(Int_t iCascSelIdx = 0; iCascSelIdx < 8; iCascSelIdx++ ) { fCascSels        [iCascSelIdx ] = -1.; }
 282 }
 283
 284
 285
 286
 287 //_____Non-default Constructor________________________________________________________________
 288 AliAnalysisTaskCheckPerformanceCascadePbPb::AliAnalysisTaskCheckPerformanceCascadePbPb(const char *name)
 289   : AliAnalysisTaskSE(name),
 290     fDebugCascade(0), fAnalysisType("ESD"), fCollidingSystems(0),
 291  fESDtrackCuts(0), /*fPaveTextBookKeeping(0),*/
 292     fkRerunV0CascVertexers         (0),
 293     fkQualityCutZprimVtxPos        (kTRUE),
 294     fkRejectEventPileUp            (kTRUE),
 295     fkQualityCutNoTPConlyPrimVtx   (kTRUE),
 296     fkQualityCutTPCrefit           (kTRUE),
 297     fkQualityCut80TPCcls           (kTRUE),
 298     fkIsDataRecoWith1PadTPCCluster (kTRUE),
 299     fkExtraSelections              (0),
 300     fUseCFCont(0),
 301     fCentrLowLim(0), fCentrUpLim(0), fCentrEstimator(0),
 302     fVtxRange                      (0),
 303     fApplyAccCut                   (0),
 304
 305
 306         // - Cascade part initialisation
 307     fListHistCascade(0),
 308
 309     // Events in centraity bins
 310     fHistEvtsInCentralityBinsvsNtracks(0),
 311
 312     // Cascade multiplicity histos
 313
 314     fHistnXiPlusPerEvTot(0),
 315     fHistnXiMinusPerEvTot(0),
 316     fHistnOmegaPlusPerEvTot(0),
 317     fHistnOmegaMinusPerEvTot(0),
 318
 319     fHistnXiPlusPerEv(0),
 320     fHistnXiMinusPerEv(0),
 321     fHistnOmegaPlusPerEv(0),
 322     fHistnOmegaMinusPerEv(0),
 323
 324     fHistnAssoXiMinus(0),
 325     fHistnAssoXiPlus(0),
 326     fHistnAssoOmegaMinus(0),
 327     fHistnAssoOmegaPlus(0),
 328
 329     fHistMCTrackMultiplicity(0),
 330        // - Resolution of the multiplicity estimator
 331     f2dHistRecoMultVsMCMult(0),
 332
 333     fHistEtaGenProton(0),
 334     fHistEtaGenAntiProton(0),
 335
 336 // Xi-
 337    fHistEtaGenCascXiMinus(0),
 338    f3dHistGenPtVsGenYGenvsCentXiMinus(0),
 339    f3dHistGenPtVsGenYGenvsNtracksXiMinus(0),
 340
 341
 342     fHistThetaGenCascXiMinus(0),
 343     f2dHistGenPtVsGenYFdblXiMinus(0),
 344
 345     fHistThetaLambdaXiMinus(0),
 346     fHistThetaBachXiMinus(0),
 347
 348     fHistThetaMesDghterXiMinus(0),
 349     fHistThetaBarDghterXiMinus(0),
 350
 351     fHistPtBachXiMinus(0),
 352     fHistPtMesDghterXiMinus(0),
 353     fHistPtBarDghterXiMinus(0),
 354
 355
 356    // Xi+
 357    fHistEtaGenCascXiPlus(0),
 358   f3dHistGenPtVsGenYGenvsCentXiPlus(0),
 359   f3dHistGenPtVsGenYGenvsNtracksXiPlus(0),
 360
 361     fHistThetaGenCascXiPlus(0),
 362     f2dHistGenPtVsGenYFdblXiPlus(0),
 363
 364     fHistThetaLambdaXiPlus(0),
 365     fHistThetaBachXiPlus(0),
 366
 367     fHistThetaMesDghterXiPlus(0),
 368     fHistThetaBarDghterXiPlus(0),
 369
 370     fHistPtBachXiPlus(0),
 371     fHistPtMesDghterXiPlus(0),
 372     fHistPtBarDghterXiPlus(0),
 373
 374    // Omega-
 375    fHistEtaGenCascOmegaMinus(0),
 376    f3dHistGenPtVsGenYGenvsCentOmegaMinus(0),
 377    f3dHistGenPtVsGenYGenvsNtracksOmegaMinus(0),
 378
 379     fHistThetaGenCascOmegaMinus(0),
 380     f2dHistGenPtVsGenYFdblOmegaMinus(0),
 381
 382     fHistThetaLambdaOmegaMinus(0),
 383     fHistThetaBachOmegaMinus(0),
 384
 385     fHistThetaMesDghterOmegaMinus(0),
 386     fHistThetaBarDghterOmegaMinus(0),
 387
 388     fHistPtBachOmegaMinus(0),
 389     fHistPtMesDghterOmegaMinus(0),
 390     fHistPtBarDghterOmegaMinus(0),
 391
 392
 393    // Omega+
 394    fHistEtaGenCascOmegaPlus(0),
 395    f3dHistGenPtVsGenYGenvsCentOmegaPlus(0),
 396    f3dHistGenPtVsGenYGenvsNtracksOmegaPlus(0),
 397
 398
 399     fHistThetaGenCascOmegaPlus(0),
 400     f2dHistGenPtVsGenYFdblOmegaPlus(0),
 401
 402     fHistThetaLambdaOmegaPlus(0),
 403     fHistThetaBachOmegaPlus(0),
 404
 405     fHistThetaMesDghterOmegaPlus(0),
 406     fHistThetaBarDghterOmegaPlus(0),
 407
 408     fHistPtBachOmegaPlus(0),
 409     fHistPtMesDghterOmegaPlus(0),
 410     fHistPtBarDghterOmegaPlus(0),
 411
 412 // Part 2 - Association to MC
 413
 414     fHistMassXiMinus(0),
 415     fHistMassXiPlus(0),
 416     fHistMassOmegaMinus(0),
 417     fHistMassOmegaPlus(0),
 418
 419         // - Effective mass histos with combined PID
 420     fHistMassWithCombPIDXiMinus(0), fHistMassWithCombPIDXiPlus(0),
 421     fHistMassWithCombPIDOmegaMinus(0), fHistMassWithCombPIDOmegaPlus(0),
 422
 423         // - PID Probability versus MC Pt(bachelor track)
 424     f2dHistPIDprobaKaonVsMCPtBach(0), f2dHistPIDprobaPionVsMCPtBach(0),
 425
 426         // - Effective mass histos with perfect MC PID on the bachelor
 427     fHistMassWithMcPIDXiMinus(0), fHistMassWithMcPIDXiPlus(0),
 428     fHistMassWithMcPIDOmegaMinus(0), fHistMassWithMcPIDOmegaPlus(0),
 429
 430
 431         // - Effective mass histos for the cascade candidates associated with MC
 432     fHistAsMCMassXiMinus(0),
 433     fHistAsMCMassXiPlus(0),
 434     fHistAsMCMassOmegaMinus(0),
 435     fHistAsMCMassOmegaPlus(0),
 436
 437         // - Generated Pt Vs generated y, for the cascade candidates associated with MC + Info Comb. PID
 438     f2dHistAsMCandCombPIDGenPtVsGenYXiMinus(0),
 439     f2dHistAsMCandCombPIDGenPtVsGenYXiPlus(0),
 440     f2dHistAsMCandCombPIDGenPtVsGenYOmegaMinus(0),
 441     f2dHistAsMCandCombPIDGenPtVsGenYOmegaPlus(0),
 442
 443         // - Generated Pt Vs generated y, for the cascade candidates associated with MC
 444     f2dHistAsMCGenPtVsGenYXiMinus(0),
 445     f2dHistAsMCGenPtVsGenYXiPlus(0),
 446     f2dHistAsMCGenPtVsGenYOmegaMinus(0),
 447     f2dHistAsMCGenPtVsGenYOmegaPlus(0),
 448
 449         // - Generated Eta of the the cascade candidates associated with MC
 450     fHistAsMCGenEtaXiMinus(0),
 451     fHistAsMCGenEtaXiPlus(0),
 452     fHistAsMCGenEtaOmegaMinus(0),
 453     fHistAsMCGenEtaOmegaPlus(0),
 454
 455         // - Resolution in Pt as function of generated Pt
 456     f2dHistAsMCResPtXiMinus(0),
 457     f2dHistAsMCResPtXiPlus(0),
 458     f2dHistAsMCResPtOmegaMinus(0),
 459     f2dHistAsMCResPtOmegaPlus(0),
 460
 461         // - Resolution in R(2D) as function of generated R
 462     f2dHistAsMCResRXiMinus(0),
 463     f2dHistAsMCResRXiPlus(0),
 464     f2dHistAsMCResROmegaMinus(0),
 465     f2dHistAsMCResROmegaPlus(0),
 466
 467        // - Resolution in phi as function of generated Pt
 468     f2dHistAsMCResPhiXiMinus(0),
 469     f2dHistAsMCResPhiXiPlus(0),
 470     f2dHistAsMCResPhiOmegaMinus(0),
 471     f2dHistAsMCResPhiOmegaPlus(0),
 472
 473     fCFContCascadePIDAsXiMinus(0),
 474     fCFContCascadePIDAsXiPlus(0),
 475     fCFContCascadePIDAsOmegaMinus(0),
 476     fCFContCascadePIDAsOmegaPlus(0),
 477
 478     fCFContAsCascadeCuts(0),
 479
 480     fV0Ampl(0)
 481
 482
 483 {
 484   // Constructor
 485
 486   // Define input and output slots here
 487   // Input slot #0 works with a TChain
 488   // Output slot #1 writes into a TList container (cascade)
 489
 490         for(Int_t iAlephIdx   = 0; iAlephIdx   < 5; iAlephIdx++   ) { fAlephParameters [iAlephIdx]    = -1.; }
 491
 492         // Hyper Loose  // FIXME change with PbPb cuts
 493
 494         fV0Sels[0] =  33.  ;  // max allowed chi2
 495         fV0Sels[1] =   0.001; // min allowed impact parameter for the 1st daughter (LHC09a4 : 0.05)
 496         fV0Sels[2] =   0.001; // min allowed impact parameter for the 2nd daughter (LHC09a4 : 0.05)
 497         fV0Sels[3] =   5.0 ;  // max allowed DCA between the daughter tracks       (LHC09a4 : 0.5)
 498         fV0Sels[4] =   0.0 ;  // min allowed cosine of V0's pointing angle         (LHC09a4 : 0.99)
 499         fV0Sels[5] =   0.1 ;  // min radius of the fiducial volume                 (LHC09a4 : 0.2)
 500         fV0Sels[6] = 100.  ;  // max radius of the fiducial volume                 (LHC09a4 : 100.0)
 501
 502         fCascSels[0] =  33.   ;  // max allowed chi2 (same as PDC07)
 503         fCascSels[1] =   0.001;  // min allowed V0 impact parameter                    (PDC07 : 0.05   / LHC09a4 : 0.025 )
 504         fCascSels[2] =   0.008;  // "window" around the Lambda mass                    (PDC07 : 0.008  / LHC09a4 : 0.010 )
 505         fCascSels[3] =   0.001;  // min allowed bachelor's impact parameter            (PDC07 : 0.035  / LHC09a4 : 0.025 )
 506         fCascSels[4] =   5.0  ;  // max allowed DCA between the V0 and the bachelor    (PDC07 : 0.1    / LHC09a4 : 0.2   )
 507         fCascSels[5] =   0.82 ;  //FIXME min allowed cosine of the cascade pointing angle   (PDC07 : 0.9985 / LHC09a4 : 0.998 )
 508         fCascSels[6] =   0.1  ;  // min radius of the fiducial volume                  (PDC07 : 0.9    / LHC09a4 : 0.2   )
 509         fCascSels[7] = 100.   ;  // max radius of the fiducial volume                  (PDC07 : 100    / LHC09a4 : 100   )
 510
 511
 512   DefineOutput(1, TList::Class());
 513
 514 }
 515
 516
 517 AliAnalysisTaskCheckPerformanceCascadePbPb::~AliAnalysisTaskCheckPerformanceCascadePbPb()
 518 {
 519   //
 520   // Destructor
 521   //
 522
 523   // For all TH1, 2, 3 HnSparse and CFContainer are in the fListCascade TList.
 524   // They will be deleted when fListCascade is deleted by the TSelector dtor
 525   // Because of TList::SetOwner()
 526
 527   if (fListHistCascade && !AliAnalysisManager::GetAnalysisManager()->IsProofMode())      { delete fListHistCascade;     fListHistCascade = 0x0;  }
 528   if (fESDtrackCuts)         { delete fESDtrackCuts;        fESDtrackCuts = 0x0; }
 529   /*if (fPaveTextBookKeeping)  { delete fPaveTextBookKeeping; fPaveTextBookKeeping = 0x0; } // fPaveTextBookKeeping is not stored into the TList*/
 530 }
 531
 532
 533 //________________________________________________________________________
 534 void AliAnalysisTaskCheckPerformanceCascadePbPb::UserCreateOutputObjects()
 535 {
 536   // Create histograms
 537   // Called once
 538
 539    // Option for AliLog
 540         AliLog::SetGlobalLogLevel(AliLog::kError);
 541         // to suppress the extensive info prompted by a run with MC
 542
 543    // Definition of the output datamembers
 544    fListHistCascade = new TList();
 545    fListHistCascade->SetOwner(); // See http://root.cern.ch/root/html/TCollection.html#TCollection:SetOwner
 546
 547    // New PID object
 548    AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 549    AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 550    fPIDResponse = inputHandler->GetPIDResponse();
 551
 552 // Only used to get the number of primary reconstructed tracks
 553 if(! fESDtrackCuts ){
 554       fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE); // Std definition of primary (see kTRUE argument) tracks for 2010
 555 //      fESDtrackCuts->SetEtaRange(-0.8,+0.8);
 556 //      fESDtrackCuts->SetPtRange(0.15, 1e10);
 557       Printf("CheckCascade - ESDtrackCuts set up to 2010 std ITS-TPC cuts...");
 558 }
 559
 560
 561 /*
 562 if( !fPaveTextBookKeeping){
 563         fPaveTextBookKeeping = new TPaveText(0.1, 0.1, 0.9, 0.9,"NDC");
 564         fPaveTextBookKeeping->SetName("fPaveTextBookKeeping");
 565         fPaveTextBookKeeping->SetBorderSize(0);
 566         fPaveTextBookKeeping->SetTextAlign(12);
 567         fPaveTextBookKeeping->SetFillColor(kWhite);
 568         fPaveTextBookKeeping->SetTextFont(42);        // regular Arial or Helvetica,
 569         fPaveTextBookKeeping->SetTextColor(kGray+3);
 570
 571
 572         fPaveTextBookKeeping->AddText( "Task CHECK PERFORMANCE CASCADE analysis" );
 573         fPaveTextBookKeeping->AddText("- - - - - - - - - - - ");
 574         fPaveTextBookKeeping->AddText( Form("AnalysisType : %s ", fAnalysisType.Data() ));
 575         if(!fCollidingSystems)  fPaveTextBookKeeping->AddText("Colliding system : p-p collisions ");
 576         else                    fPaveTextBookKeeping->AddText("Colliding system : A-A collisions ");
 577
 578         fPaveTextBookKeeping->AddText("- - - - - - - - - - - ");
 579
 580         if(fkRerunV0CascVertexers){
 581                 fPaveTextBookKeeping->AddText("A.1. With V0 vertexer : ");
 582                 fPaveTextBookKeeping->AddText( Form("  - V0 #chi^{2} _________________ <  %.3f ",               fV0Sels[0] ));
 583                 fPaveTextBookKeeping->AddText( Form("  - DCA(prim. Vtx/ 1^{st} daughter) ___ >  %.3f     cm ",  fV0Sels[1] ));
 584                 fPaveTextBookKeeping->AddText( Form("  - DCA(prim. Vtx/ 2^{nd} daughter) __  >  %.3f     cm",   fV0Sels[2] ));
 585                 fPaveTextBookKeeping->AddText( Form("  - DCA between V0 daughters ___ <  %.3f      cm",         fV0Sels[3] ));
 586                 fPaveTextBookKeeping->AddText( Form("  - cos(V0 pointing angle) _______ >  %.3f ",              fV0Sels[4] ));
 587                 fPaveTextBookKeeping->AddText( Form("  - R_{transv}(V0 decay) ________ >  %.3f             cm", fV0Sels[5] ));
 588                 fPaveTextBookKeeping->AddText( Form("  - R_{transv}(V0 decay) ________ <  %.3f         cm",     fV0Sels[6] ));
 589
 590                 fPaveTextBookKeeping->AddText(" ");
 591
 592                 fPaveTextBookKeeping->AddText("A.2. With Casc. vertexer : ");
 593                 fPaveTextBookKeeping->AddText( Form("  - Casc. #chi^{2} ______________  <  %.3f ",                               fCascSels[0] ));
 594                 fPaveTextBookKeeping->AddText( Form("  - DCA(prim. Vtx/ V0) _________ >  %.3f    cm",                            fCascSels[1] ));
 595                 fPaveTextBookKeeping->AddText( Form("  - | M_{#Lambda}(reco) - M_{#Lambda}(pdg) | _______ <  %.3f    GeV/c^{2}", fCascSels[2] ));
 596                 fPaveTextBookKeeping->AddText( Form("  - DCA(prim. Vtx/ Bach) _______ >  %.3f    cm",                            fCascSels[3] ));
 597                 fPaveTextBookKeeping->AddText( Form("  - DCA between Bach/ #Lambda ______ <  %.3f    cm",                        fCascSels[4] ));
 598                 fPaveTextBookKeeping->AddText( Form("  - cos(Casc. pointing angle) ____ >  %.3f ",                               fCascSels[5] ));
 599                 fPaveTextBookKeeping->AddText( Form("  - R_{transv}(Casc. decay) ______ >  %.3f         cm",                     fCascSels[6] ));
 600                 fPaveTextBookKeeping->AddText( Form("  - R_{transv}(Casc. decay) ______ <  %.3f     cm",                         fCascSels[7] ));
 601         }
 602         else{   fPaveTextBookKeeping->AddText("A. No rerunning of the V0/Casc. vertexers ... See std cuts in (AliRoot+Rec.C) used for this prod. cycle");}
 603
 604         fPaveTextBookKeeping->AddText("- - - - - - - - - - - ");
 605
 606         if(fkQualityCutZprimVtxPos)      fPaveTextBookKeeping->AddText("B. Quality Cut(prim. Vtx z-Pos)    = ON  ");
 607         else                             fPaveTextBookKeeping->AddText("B. Quality Cut(prim. Vtx z-Pos)    = Off ");
 608         if(fkQualityCutNoTPConlyPrimVtx) fPaveTextBookKeeping->AddText("C. Quality Cut(No TPC prim. vtx) = ON  ");
 609         else                             fPaveTextBookKeeping->AddText("C. Quality Cut(No TPC prim. vtx) = Off ");
 610         if(fkQualityCutTPCrefit)         fPaveTextBookKeeping->AddText("D. Quality Cut(TPCrefit)               = ON  ");
 611         else                             fPaveTextBookKeeping->AddText("D. Quality Cut(TPCrefit)               = Off ");
 612         if(fkQualityCut80TPCcls)         fPaveTextBookKeeping->AddText("E. Quality Cut(80 TPC clusters)   = ON  ");
 613         else                             fPaveTextBookKeeping->AddText("E. Quality Cut(80 TPC clusters)   = Off ");
 614         if(fkExtraSelections)            fPaveTextBookKeeping->AddText("F. Extra Analysis Selections         = ON  ");
 615         else                             fPaveTextBookKeeping->AddText("F. Extra Analysis Selections         = Off ");
 616
 617         fPaveTextBookKeeping->AddText("- - - - - - - - - - - ");
 618
 619         fPaveTextBookKeeping->AddText("G. TPC Aleph Param : ");
 620         fPaveTextBookKeeping->AddText( Form("   - fAlephParam [0] =  %.5g", fAlephParameters[0] ));
 621         fPaveTextBookKeeping->AddText( Form("   - fAlephParam [1] =  %.5g", fAlephParameters[1] ));
 622         fPaveTextBookKeeping->AddText( Form("   - fAlephParam [2] =  %.5g", fAlephParameters[2] ));
 623         fPaveTextBookKeeping->AddText( Form("   - fAlephParam [3] =  %.5g", fAlephParameters[3] ));
 624         fPaveTextBookKeeping->AddText( Form("   - fAlephParam [4] =  %.5g", fAlephParameters[4] ));
 625
 626         fListHistCascade->Add(fPaveTextBookKeeping);
 627 }
 628 */
 629
 630   // - General
 631   Double_t ptBinLimits[101];
 632   for (Int_t iptbin = 0; iptbin<101; ++iptbin) {ptBinLimits[iptbin]=iptbin*0.1;};
 633   Double_t yBinLimits[221];
 634   for (Int_t iybin = 0; iybin<221; ++iybin) {yBinLimits[iybin]=-1.1+iybin*0.01;};
 635
 636   // Events in centrality bins
 637   Double_t centBinLimits[12] = {0.,5.,10.,20.,30.,40.,50.,60.,70.,80.,90.,100.};
 638   fHistEvtsInCentralityBinsvsNtracks = new TH2F("fHistEvtsInCentralityBinsvsNtracks","",11,centBinLimits,100,0.,6000.);
 639   fListHistCascade->Add(fHistEvtsInCentralityBinsvsNtracks);
 640
 641   // Cascade multiplicity distributions
 642
 643   fHistnXiPlusPerEvTot= new TH1F("fHistnXiPlusPerEvTot", "", 100, 0, 100);
 644   fListHistCascade->Add(fHistnXiPlusPerEvTot);
 645   fHistnXiMinusPerEvTot= new TH1F("fHistnXiMinusPerEvTot", "", 100, 0, 100);
 646   fListHistCascade->Add(fHistnXiMinusPerEvTot);
 647   fHistnOmegaPlusPerEvTot = new TH1F("fHistnOmegaPlusPerEvTot", "", 50, 0, 50);
 648   fListHistCascade->Add(fHistnOmegaPlusPerEvTot);
 649   fHistnOmegaMinusPerEvTot= new TH1F("fHistnOmegaMinusPerEvTot", "", 50, 0, 50);
 650   fListHistCascade->Add(fHistnOmegaMinusPerEvTot);
 651
 652   fHistnXiPlusPerEv= new TH1F("fHistnXiPlusPerEv", "", 100, 0, 100);
 653   fListHistCascade->Add(fHistnXiPlusPerEv);
 654   fHistnXiMinusPerEv= new TH1F("fHistnXiMinusPerEv", "", 100, 0, 100);
 655   fListHistCascade->Add(fHistnXiMinusPerEv);
 656   fHistnOmegaPlusPerEv= new TH1F("fHistnOmegaPlusPerEv", "", 50, 0, 50);
 657   fListHistCascade->Add(fHistnOmegaPlusPerEv);
 658   fHistnOmegaMinusPerEv= new TH1F("fHistnOmegaMinusPerEv", "", 50, 0, 50);
 659   fListHistCascade->Add(fHistnOmegaMinusPerEv);
 660
 661   fHistnAssoXiMinus= new TH1F("fHistnAssoXiMinus", "", 100, 0, 100);
 662   fListHistCascade->Add(fHistnAssoXiMinus);
 663   fHistnAssoXiPlus= new TH1F("fHistnAssoXiPlus", "", 100, 0, 100);
 664   fListHistCascade->Add(fHistnAssoXiPlus);
 665   fHistnAssoOmegaMinus= new TH1F("fHistnAssoOmegaMinus", "", 50, 0, 50);
 666   fListHistCascade->Add(fHistnAssoOmegaMinus);
 667   fHistnAssoOmegaPlus= new TH1F("fHistnAssoOmegaPlus", "", 50, 0, 50);
 668   fListHistCascade->Add(fHistnAssoOmegaPlus);
 669
 670   if (!fHistMCTrackMultiplicity) {
 671      fHistMCTrackMultiplicity = new TH1F("fHistMCTrackMultiplicity", "Multiplicity distribution;Number of tracks;Events", 200, 0, 400000);
 672   fListHistCascade->Add(fHistMCTrackMultiplicity);
 673   }
 674
 675     // - Resolution of the multiplicity estimator
 676   if(! f2dHistRecoMultVsMCMult){
 677        f2dHistRecoMultVsMCMult = new TH2F("f2dHistRecoMultVsMCMult", "Resolution of the multiplicity estimator (prim. tracks); Reco Multiplicity (prim. tracks); MC multiplicity (gen. part.)", 200, 0., 6000., 200, 0., 6000.);
 678        fListHistCascade->Add(f2dHistRecoMultVsMCMult);
 679   }
 680
 681   if (!fHistEtaGenProton) {
 682      fHistEtaGenProton = new TH1F("fHistEtaGenProton", "#eta of any gen. p^{+};#eta;Number of prim. protons", 200, -10, 10);
 683      fListHistCascade->Add(fHistEtaGenProton);
 684   }
 685
 686   if (!fHistEtaGenAntiProton) {
 687      fHistEtaGenAntiProton = new TH1F("fHistEtaGenAntiProton", "#eta of any gen. #bar{p}^{-};#eta;Number of prim. #bar{p}", 200, -10, 10);
 688      fListHistCascade->Add(fHistEtaGenAntiProton);
 689   }
 690
 691
 692   //--------
 693   // I - Xi-
 694   // - Pseudo-Rapidity distribution
 695   if (!fHistEtaGenCascXiMinus) {
 696      fHistEtaGenCascXiMinus = new TH1F("fHistEtaGenCascXiMinus", "#eta of any gen. #Xi^{-};#eta;Number of Casc", 200, -10, 10);
 697      fListHistCascade->Add(fHistEtaGenCascXiMinus);
 698   }
 699
 700   if (!f3dHistGenPtVsGenYGenvsCentXiMinus) {
 701      f3dHistGenPtVsGenYGenvsCentXiMinus = new TH3D("f3dHistGenPtVsGenYGenvsCentXiMinus", "MC P_{t} Vs MC Y of Gen #Xi^{-} ;Pt_{MC} (GeV/c); Y_{MC}", 100, ptBinLimits, 220, yBinLimits, 11, centBinLimits);
 702      fListHistCascade->Add(f3dHistGenPtVsGenYGenvsCentXiMinus);
 703   }
 704   if (!f3dHistGenPtVsGenYGenvsNtracksXiMinus) {
 705      f3dHistGenPtVsGenYGenvsNtracksXiMinus = new TH3D("f3dHistGenPtVsGenYGenvsNtracksXiMinus", "MC P_{t} Vs MC Y of Gen #Xi^{-} ;Pt_{MC} (GeV/c); Y_{MC}", 100, 0., 10., 220, -1.1, 1.1, 100, 0., 6000.);
 706      fListHistCascade->Add(f3dHistGenPtVsGenYGenvsNtracksXiMinus);
 707   }
 708
 709
 710
 711   // - Info at the generation level of multi-strange particle
 712
 713   if (!fHistThetaGenCascXiMinus) {
 714      fHistThetaGenCascXiMinus = new TH1F("fHistThetaGenCascXiMinus", "#theta of gen. #Xi^{-};#theta;Number of Casc.", 200, -10, 190);
 715      fListHistCascade->Add(fHistThetaGenCascXiMinus);
 716   }
 717
 718   if (!f2dHistGenPtVsGenYFdblXiMinus) {
 719      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);
 720      fListHistCascade->Add(f2dHistGenPtVsGenYFdblXiMinus);
 721   }
 722
 723                 // - Theta distribution the daughters (control plots)
 724
 725   if (!fHistThetaLambdaXiMinus) {
 726      fHistThetaLambdaXiMinus = new TH1F("fHistThetaLambdaXiMinus", "#theta of gen. #Lambda (Xi dghter);#theta_{#Lambda};Number of #Lambda^0", 200, -10, 190);
 727      fListHistCascade->Add(fHistThetaLambdaXiMinus);
 728   }
 729
 730   if (!fHistThetaBachXiMinus) {
 731      fHistThetaBachXiMinus = new TH1F("fHistThetaBachXiMinus", "#theta of gen. Bach.;#theta_{Bach};Number of Bach.", 200, -10, 190);
 732      fListHistCascade->Add(fHistThetaBachXiMinus);
 733   }
 734
 735   if (!fHistThetaMesDghterXiMinus) {
 736      fHistThetaMesDghterXiMinus = new TH1F("fHistThetaMesDghterXiMinus", "#theta of gen. Meson #Lambda dghter;#theta_{MesDght};Number of Mes.", 200, -10, 190);
 737      fListHistCascade->Add(fHistThetaMesDghterXiMinus);
 738   }
 739
 740   if (!fHistThetaBarDghterXiMinus) {
 741      fHistThetaBarDghterXiMinus = new TH1F("fHistThetaBarDghterXiMinus", "#theta of gen. Baryon #Lambda dghter;#theta_{BarDght};Number of Bar.", 200, -10, 190);
 742      fListHistCascade->Add(fHistThetaBarDghterXiMinus);
 743   }
 744
 745                 // - Pt distribution (control plots)
 746
 747   if (!fHistPtBachXiMinus) {
 748      fHistPtBachXiMinus = new TH1F("fHistPtBachXiMinus", "p_{t} of gen. Bach.;pt_{Bach};Number of Bach.", 200, 0, 10);
 749      fListHistCascade->Add(fHistPtBachXiMinus);
 750   }
 751
 752   if (!fHistPtMesDghterXiMinus) {
 753      fHistPtMesDghterXiMinus = new TH1F("fHistPtMesDghterXiMinus", "p_{t} of gen. Meson #Lambda dghter;pt_{MesDght};Number of Mes.", 200, 0, 10);
 754      fListHistCascade->Add(fHistPtMesDghterXiMinus);
 755   }
 756
 757   if (!fHistPtBarDghterXiMinus) {
 758      fHistPtBarDghterXiMinus = new TH1F("fHistPtBarDghterXiMinus", "p_{t} of gen. Baryon #Lambda dghter;pt_{BarDght};Number of Bar.", 200, 0, 10);
 759      fListHistCascade->Add(fHistPtBarDghterXiMinus);
 760   }
 761
 762
 763
 764   //--------
 765   // II - Xi+
 766   // - Pseudo-Rapidity distribution
 767   if (!fHistEtaGenCascXiPlus) {
 768      fHistEtaGenCascXiPlus = new TH1F("fHistEtaGenCascXiPlus", "#eta of any gen. #Xi^{+};#eta;Number of Casc", 200, -10, 10);
 769      fListHistCascade->Add(fHistEtaGenCascXiPlus);
 770   }
 771 if (!f3dHistGenPtVsGenYGenvsCentXiPlus) {
 772      f3dHistGenPtVsGenYGenvsCentXiPlus = new TH3D("f3dHistGenPtVsGenYGenvsCentXiPlus", "MC P_{t} Vs MC Y of Gen #Xi^{-} ;Pt_{MC} (GeV/c); Y_{MC}", 100, ptBinLimits, 220, yBinLimits, 11, centBinLimits);
 773      fListHistCascade->Add(f3dHistGenPtVsGenYGenvsCentXiPlus);
 774   }
 775   if (!f3dHistGenPtVsGenYGenvsNtracksXiPlus) {
 776      f3dHistGenPtVsGenYGenvsNtracksXiPlus = new TH3D("f3dHistGenPtVsGenYGenvsNtracksXiPlus", "MC P_{t} Vs MC Y of Gen #Xi^{-} ;Pt_{MC} (GeV/c); Y_{MC}", 100, 0., 10., 220, -1.1, 1.1, 100, 0., 6000.);
 777      fListHistCascade->Add(f3dHistGenPtVsGenYGenvsNtracksXiPlus);
 778   }
 779
 780
 781   // - Info at the generation level of multi-strange particle
 782
 783   if (!fHistThetaGenCascXiPlus) {
 784      fHistThetaGenCascXiPlus = new TH1F("fHistThetaGenCascXiPlus", "#theta of gen. #Xi^{+};#theta;Number of Casc.", 200, -10, 190);
 785      fListHistCascade->Add(fHistThetaGenCascXiPlus);
 786   }
 787
 788   if (!f2dHistGenPtVsGenYFdblXiPlus) {
 789      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);
 790      fListHistCascade->Add(f2dHistGenPtVsGenYFdblXiPlus);
 791   }
 792
 793                 // - Theta distribution the daughters (control plots)
 794
 795   if (!fHistThetaLambdaXiPlus) {
 796      fHistThetaLambdaXiPlus = new TH1F("fHistThetaLambdaXiPlus", "#theta of gen. #Lambda (Xi dghter);#theta_{#Lambda};Number of #Lambda", 200, -10, 190);
 797      fListHistCascade->Add(fHistThetaLambdaXiPlus);
 798   }
 799
 800   if (!fHistThetaBachXiPlus) {
 801      fHistThetaBachXiPlus = new TH1F("fHistThetaBachXiPlus", "#theta of gen. Bach.;#theta_{Bach};Number of Bach.", 200, -10, 190);
 802      fListHistCascade->Add(fHistThetaBachXiPlus);
 803   }
 804
 805   if (!fHistThetaMesDghterXiPlus) {
 806      fHistThetaMesDghterXiPlus = new TH1F("fHistThetaMesDghterXiPlus", "#theta of gen. Meson #Lambda dghter;#theta_{MesDght};Number of Mes.", 200, -10, 190);
 807      fListHistCascade->Add(fHistThetaMesDghterXiPlus);
 808   }
 809
 810   if (!fHistThetaBarDghterXiPlus) {
 811      fHistThetaBarDghterXiPlus = new TH1F("fHistThetaBarDghterXiPlus", "#theta of gen. Baryon #Lambda dghter;#theta_{BarDght};Number of Bar.", 200, -10, 190);
 812      fListHistCascade->Add(fHistThetaBarDghterXiPlus);
 813   }
 814
 815                 // - Pt distribution (control plots)
 816
 817   if (!fHistPtBachXiPlus) {
 818      fHistPtBachXiPlus = new TH1F("fHistPtBachXiPlus", "p_{t} of gen. Bach.;pt_{Bach};Number of Bach.", 200, 0, 10);
 819      fListHistCascade->Add(fHistPtBachXiPlus);
 820   }
 821
 822   if (!fHistPtMesDghterXiPlus) {
 823      fHistPtMesDghterXiPlus = new TH1F("fHistPtMesDghterXiPlus", "p_{t} of gen. Meson #Lambda dghter);pt_{MesDght};Number of Mes.", 200, 0, 10);
 824      fListHistCascade->Add(fHistPtMesDghterXiPlus);
 825   }
 826
 827   if (!fHistPtBarDghterXiPlus) {
 828      fHistPtBarDghterXiPlus = new TH1F("fHistPtBarDghterXiPlus", "p_{t} of gen. Baryon #Lambda dghter);pt_{BarDght};Number of Bar.", 200, 0, 10);
 829      fListHistCascade->Add(fHistPtBarDghterXiPlus);
 830   }
 831
 832
 833   //---------
 834   // III - Omega-
 835   // - Pseudo-Rapidity distribution
 836   if (!fHistEtaGenCascOmegaMinus) {
 837      fHistEtaGenCascOmegaMinus = new TH1F("fHistEtaGenCascOmegaMinus", "#eta of any gen. #Omega^{-};#eta;Number of Casc", 200, -10, 10);
 838      fListHistCascade->Add(fHistEtaGenCascOmegaMinus);
 839   }
 840  if (!f3dHistGenPtVsGenYGenvsCentOmegaMinus) {
 841      f3dHistGenPtVsGenYGenvsCentOmegaMinus = new TH3D("f3dHistGenPtVsGenYGenvsCentOmegaMinus", "MC P_{t} Vs MC Y of Gen #Xi^{-} ;Pt_{MC} (GeV/c); Y_{MC}", 100, ptBinLimits, 220, yBinLimits, 11, centBinLimits);
 842      fListHistCascade->Add(f3dHistGenPtVsGenYGenvsCentOmegaMinus);
 843   }
 844   if (!f3dHistGenPtVsGenYGenvsNtracksOmegaMinus) {
 845      f3dHistGenPtVsGenYGenvsNtracksOmegaMinus = new TH3D("f3dHistGenPtVsGenYGenvsNtracksOmegaMinus", "MC P_{t} Vs MC Y of Gen #Xi^{-} ;Pt_{MC} (GeV/c); Y_{MC}", 100, 0., 10., 220, -1.1, 1.1, 100, 0., 6000.);
 846      fListHistCascade->Add(f3dHistGenPtVsGenYGenvsNtracksOmegaMinus);
 847   }
 848
 849
 850
 851   // - Info at the generation level of multi-strange particle
 852
 853   if (!fHistThetaGenCascOmegaMinus) {
 854      fHistThetaGenCascOmegaMinus = new TH1F("fHistThetaGenCascOmegaMinus", "#theta of gen. #Omega^{-};#theta;Number of Casc.", 200, -10, 190);
 855      fListHistCascade->Add(fHistThetaGenCascOmegaMinus);
 856   }
 857
 858   if (!f2dHistGenPtVsGenYFdblOmegaMinus) {
 859      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);
 860      fListHistCascade->Add(f2dHistGenPtVsGenYFdblOmegaMinus);
 861   }
 862
 863                 // - Theta distribution the daughters (control plots)
 864
 865   if (!fHistThetaLambdaOmegaMinus) {
 866      fHistThetaLambdaOmegaMinus = new TH1F("fHistThetaLambdaOmegaMinus", "#theta of gen. #Lambda (Omega dghter);#theta_{#Lambda};Number of #Lambda", 200, -10, 190);
 867      fListHistCascade->Add(fHistThetaLambdaOmegaMinus);
 868   }
 869
 870   if (!fHistThetaBachOmegaMinus) {
 871      fHistThetaBachOmegaMinus = new TH1F("fHistThetaBachOmegaMinus", "#theta of gen. Bach.;#theta_{Bach};Number of Bach.", 200, -10, 190);
 872      fListHistCascade->Add(fHistThetaBachOmegaMinus);
 873   }
 874
 875   if (!fHistThetaMesDghterOmegaMinus) {
 876      fHistThetaMesDghterOmegaMinus = new TH1F("fHistThetaMesDghterOmegaMinus", "#theta of gen. Meson #Lambda dghter;#theta_{MesDght};Number of Mes.", 200, -10, 190);
 877      fListHistCascade->Add(fHistThetaMesDghterOmegaMinus);
 878   }
 879
 880   if (!fHistThetaBarDghterOmegaMinus) {
 881      fHistThetaBarDghterOmegaMinus = new TH1F("fHistThetaBarDghterOmegaMinus", "#theta of gen. Baryon #Lambda dghter;#theta_{BarDght};Number of Bar.", 200, -10, 190);
 882      fListHistCascade->Add(fHistThetaBarDghterOmegaMinus);
 883   }
 884
 885                 // - Pt distribution (control plots)
 886
 887   if (!fHistPtBachOmegaMinus) {
 888      fHistPtBachOmegaMinus = new TH1F("fHistPtBachOmegaMinus", "p_{t} of gen. Bach.;pt_{Bach};Number of Bach.", 200, 0, 10);
 889      fListHistCascade->Add(fHistPtBachOmegaMinus);
 890   }
 891
 892   if (!fHistPtMesDghterOmegaMinus) {
 893      fHistPtMesDghterOmegaMinus = new TH1F("fHistPtMesDghterOmegaMinus", "p_{t} of gen. Meson #Lambda dghter);pt_{MesDght};Number of Mes.", 200, 0, 10);
 894      fListHistCascade->Add(fHistPtMesDghterOmegaMinus);
 895   }
 896
 897   if (!fHistPtBarDghterOmegaMinus) {
 898      fHistPtBarDghterOmegaMinus = new TH1F("fHistPtBarDghterOmegaMinus", "p_{t} of gen. Baryon #Lambda dghter);pt_{BarDght};Number of Bar.", 200, 0, 10);
 899      fListHistCascade->Add(fHistPtBarDghterOmegaMinus);
 900   }
 901
 902
 903   //---------
 904   // IV - Omega+
 905   // - Pseudo-Rapidity distribution
 906   if (!fHistEtaGenCascOmegaPlus) {
 907      fHistEtaGenCascOmegaPlus = new TH1F("fHistEtaGenCascOmegaPlus", "#eta of any gen. #Omega^{+};#eta;Number of Casc", 200, -10, 10);
 908      fListHistCascade->Add(fHistEtaGenCascOmegaPlus);
 909   }
 910   if (!f3dHistGenPtVsGenYGenvsCentOmegaPlus) {
 911      f3dHistGenPtVsGenYGenvsCentOmegaPlus = new TH3D("f3dHistGenPtVsGenYGenvsCentOmegaPlus", "MC P_{t} Vs MC Y of Gen #Xi^{-} ;Pt_{MC} (GeV/c); Y_{MC}", 100, ptBinLimits, 220, yBinLimits, 11, centBinLimits);
 912      fListHistCascade->Add(f3dHistGenPtVsGenYGenvsCentOmegaPlus);
 913   }
 914   if (!f3dHistGenPtVsGenYGenvsNtracksOmegaPlus) {
 915      f3dHistGenPtVsGenYGenvsNtracksOmegaPlus = new TH3D("f3dHistGenPtVsGenYGenvsNtracksOmegaPlus", "MC P_{t} Vs MC Y of Gen #Xi^{-} ;Pt_{MC} (GeV/c); Y_{MC}", 100, 0., 10., 220, -1.1, 1.1, 100, 0., 6000.);
 916      fListHistCascade->Add(f3dHistGenPtVsGenYGenvsNtracksOmegaPlus);
 917   }
 918
 919
 920   // - Info at the generation level of multi-strange particle
 921
 922   if (!fHistThetaGenCascOmegaPlus) {
 923      fHistThetaGenCascOmegaPlus = new TH1F("fHistThetaGenCascOmegaPlus", "#theta of gen. #Omega^{+};#theta;Number of Casc.", 200, -10, 190);
 924      fListHistCascade->Add(fHistThetaGenCascOmegaPlus);
 925   }
 926
 927   if (!f2dHistGenPtVsGenYFdblOmegaPlus) {
 928      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);
 929      fListHistCascade->Add(f2dHistGenPtVsGenYFdblOmegaPlus);
 930   }
 931
 932
 933                 // - Theta distribution the daughters (control plots)
 934
 935   if (!fHistThetaLambdaOmegaPlus) {
 936      fHistThetaLambdaOmegaPlus = new TH1F("fHistThetaLambdaOmegaPlus", "#theta of gen. #Lambda (Omega dghter);#theta_{#Lambda};Number of #Lambda", 200, -10, 190);
 937      fListHistCascade->Add(fHistThetaLambdaOmegaPlus);
 938   }
 939
 940   if (!fHistThetaBachOmegaPlus) {
 941      fHistThetaBachOmegaPlus = new TH1F("fHistThetaBachOmegaPlus", "#theta of gen. Bach.;#theta_{Bach};Number of Bach.", 200, -10, 190);
 942      fListHistCascade->Add(fHistThetaBachOmegaPlus);
 943   }
 944
 945   if (!fHistThetaMesDghterOmegaPlus) {
 946      fHistThetaMesDghterOmegaPlus = new TH1F("fHistThetaMesDghterOmegaPlus", "#theta of gen. Meson #Lambda dghter;#theta_{MesDght};Number of Mes.", 200, -10, 190);
 947      fListHistCascade->Add(fHistThetaMesDghterOmegaPlus);
 948   }
 949
 950   if (!fHistThetaBarDghterOmegaPlus) {
 951      fHistThetaBarDghterOmegaPlus = new TH1F("fHistThetaBarDghterOmegaPlus", "#theta of gen. Baryon #Lambda dghter;#theta_{BarDght};Number of Bar.", 200, -10, 190);
 952      fListHistCascade->Add(fHistThetaBarDghterOmegaPlus);
 953   }
 954
 955                 // - Pt distribution (control plots)
 956
 957   if (!fHistPtBachOmegaPlus) {
 958      fHistPtBachOmegaPlus = new TH1F("fHistPtBachOmegaPlus", "p_{t} of gen. Bach.;pt_{Bach};Number of Bach.", 200, 0, 10);
 959      fListHistCascade->Add(fHistPtBachOmegaPlus);
 960   }
 961
 962   if (!fHistPtMesDghterOmegaPlus) {
 963      fHistPtMesDghterOmegaPlus = new TH1F("fHistPtMesDghterOmegaPlus", "p_{t} of gen. Meson #Lambda dghter);pt_{MesDght};Number of Mes.", 200, 0, 10);
 964      fListHistCascade->Add(fHistPtMesDghterOmegaPlus);
 965   }
 966
 967   if (!fHistPtBarDghterOmegaPlus) {
 968      fHistPtBarDghterOmegaPlus = new TH1F("fHistPtBarDghterOmegaPlus", "p_{t} of gen. Baryon #Lambda dghter);pt_{BarDght};Number of Bar.", 200, 0, 10);
 969      fListHistCascade->Add(fHistPtBarDghterOmegaPlus);
 970   }
 971
 972
 973 //--------------------------------------------------------------------------------
 974 // Part 2 - Any reconstructed cascades + reconstructed cascades associated with MC
 975
 976                 // - Effective mass histos for cascades candidates.
 977
 978   if (! fHistMassXiMinus) {
 979           fHistMassXiMinus = new TH1F("fHistMassXiMinus","#Xi^{-} candidates;M( #Lambda , #pi^{-} ) (GeV/c^{2});Counts", 400,1.2,2.0);
 980           fListHistCascade->Add(fHistMassXiMinus);
 981   }
 982
 983   if (! fHistMassXiPlus) {
 984           fHistMassXiPlus = new TH1F("fHistMassXiPlus","#Xi^{+} candidates;M( #bar{#Lambda}^{0} , #pi^{+} ) (GeV/c^{2});Counts",400,1.2,2.0);
 985           fListHistCascade->Add(fHistMassXiPlus);
 986   }
 987
 988   if (! fHistMassOmegaMinus) {
 989           fHistMassOmegaMinus = new TH1F("fHistMassOmegaMinus","#Omega^{-} candidates;M( #Lambda , K^{-} ) (GeV/c^{2});Counts", 500,1.5,2.5);
 990           fListHistCascade->Add(fHistMassOmegaMinus);
 991   }
 992
 993   if (! fHistMassOmegaPlus) {
 994           fHistMassOmegaPlus = new TH1F("fHistMassOmegaPlus","#Omega^{+} candidates;M( #bar{#Lambda}^{0} , K^{+} ) (GeV/c^{2});Counts", 500,1.5,2.5);
 995           fListHistCascade->Add(fHistMassOmegaPlus);
 996   }
 997
 998
 999
1000                 // - Effective mass histos with combined PID
1001
1002   if (! fHistMassWithCombPIDXiMinus) {
1003     fHistMassWithCombPIDXiMinus = new TH1F("fHistMassWithCombPIDXiMinus","#Xi^{-} candidates, with Bach. comb. PID;M( #Lambda , #pi^{-} ) (GeV/c^{2});Counts", 400,1.2,2.0);
1004     fListHistCascade->Add(fHistMassWithCombPIDXiMinus);
1005   }
1006
1007   if (! fHistMassWithCombPIDXiPlus) {
1008     fHistMassWithCombPIDXiPlus = new TH1F("fHistMassWithCombPIDXiPlus","#Xi^{+} candidates, with Bach. comb. PID;M( #bar{#Lambda}^{0} , #pi^{+} ) (GeV/c^{2});Counts",400,1.2,2.0);
1009     fListHistCascade->Add(fHistMassWithCombPIDXiPlus);
1010   }
1011
1012   if (! fHistMassWithCombPIDOmegaMinus) {
1013         fHistMassWithCombPIDOmegaMinus = new TH1F("fHistMassWithCombPIDOmegaMinus","#Omega^{-} candidates, with Bach. comb. PID;M( #Lambda , K^{-} ) (GeV/c^{2});Counts", 500,1.5,2.5);
1014     fListHistCascade->Add(fHistMassWithCombPIDOmegaMinus);
1015   }
1016
1017   if (! fHistMassWithCombPIDOmegaPlus) {
1018         fHistMassWithCombPIDOmegaPlus = new TH1F("fHistMassWithCombPIDOmegaPlus","#Omega^{+} candidates, with Bach. comb. PID;M( #bar{#Lambda}^{0} , K^{+} ) (GeV/c^{2});Counts", 500,1.5,2.5);
1019     fListHistCascade->Add(fHistMassWithCombPIDOmegaPlus);
1020   }
1021
1022                 // - PID Probability versus MC Pt(bachelor track)
1023   if(! f2dHistPIDprobaKaonVsMCPtBach ){
1024         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 );
1025         fListHistCascade->Add(f2dHistPIDprobaKaonVsMCPtBach);
1026   }
1027
1028   if(! f2dHistPIDprobaPionVsMCPtBach ){
1029         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 );
1030         fListHistCascade->Add(f2dHistPIDprobaPionVsMCPtBach);
1031   }
1032
1033
1034                 // - Effective mass histos with perfect MC PID on the bachelor
1035
1036   if (! fHistMassWithMcPIDXiMinus) {
1037     fHistMassWithMcPIDXiMinus = new TH1F("fHistMassWithMcPIDXiMinus","#Xi^{-} candidates, with Bach. MC PID;M( #Lambda , #pi^{-} ) (GeV/c^{2});Counts", 400,1.2,2.0);
1038     fListHistCascade->Add(fHistMassWithMcPIDXiMinus);
1039   }
1040
1041   if (! fHistMassWithMcPIDXiPlus) {
1042     fHistMassWithMcPIDXiPlus = new TH1F("fHistMassWithMcPIDXiPlus","#Xi^{+} candidates, with Bach. MC PID;M( #bar{#Lambda}^{0} , #pi^{+} ) (GeV/c^{2});Counts",400,1.2,2.0);
1043     fListHistCascade->Add(fHistMassWithMcPIDXiPlus);
1044   }
1045
1046   if (! fHistMassWithMcPIDOmegaMinus) {
1047         fHistMassWithMcPIDOmegaMinus = new TH1F("fHistMassWithMcPIDOmegaMinus","#Omega^{-} candidates, with Bach. MC PID;M( #Lambda , K^{-} ) (GeV/c^{2});Counts", 500,1.5,2.5);
1048     fListHistCascade->Add(fHistMassWithMcPIDOmegaMinus);
1049   }
1050
1051   if (! fHistMassWithMcPIDOmegaPlus) {
1052         fHistMassWithMcPIDOmegaPlus = new TH1F("fHistMassWithMcPIDOmegaPlus","#Omega^{+} candidates, with Bach. MC PID;M( #bar{#Lambda}^{0} , K^{+} ) (GeV/c^{2});Counts", 500,1.5,2.5);
1053     fListHistCascade->Add(fHistMassWithMcPIDOmegaPlus);
1054   }
1055
1056
1057                 // - Effective mass histos for cascades candidates ASSOCIATED with MC.
1058
1059   if (! fHistAsMCMassXiMinus) {
1060           fHistAsMCMassXiMinus = new TH1F("fHistAsMCMassXiMinus","#Xi^{-} candidates associated to MC;M( #Lambda , #pi^{-} ) (GeV/c^{2});Counts", 400,1.2,2.0);
1061           fListHistCascade->Add(fHistAsMCMassXiMinus);
1062   }
1063
1064   if (! fHistAsMCMassXiPlus) {
1065           fHistAsMCMassXiPlus = new TH1F("fHistAsMCMassXiPlus","#Xi^{+} candidates associated to MC;M( #bar{#Lambda}^{0} , #pi^{+} ) (GeV/c^{2});Counts",400,1.2,2.0);
1066           fListHistCascade->Add(fHistAsMCMassXiPlus);
1067   }
1068
1069   if (! fHistAsMCMassOmegaMinus) {
1070           fHistAsMCMassOmegaMinus = new TH1F("fHistAsMCMassOmegaMinus","#Omega^{-} candidates associated to MC;M( #Lambda , K^{-} ) (GeV/c^{2});Counts", 500,1.5,2.5);
1071           fListHistCascade->Add(fHistAsMCMassOmegaMinus);
1072   }
1073
1074   if (! fHistAsMCMassOmegaPlus) {
1075           fHistAsMCMassOmegaPlus = new TH1F("fHistAsMCMassOmegaPlus","#Omega^{+} candidates associated to MC;M( #bar{#Lambda}^{0} , K^{+} ) (GeV/c^{2});Counts", 500,1.5,2.5);
1076           fListHistCascade->Add(fHistAsMCMassOmegaPlus);
1077   }
1078
1079
1080                 // -  Generated Pt Vs generated Y of the cascade candidates associated with MC
1081                 //     + having the proper maximum proba of combined PID for the bachelor
1082
1083   if (!f2dHistAsMCandCombPIDGenPtVsGenYXiMinus) {
1084      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);
1085      fListHistCascade->Add(f2dHistAsMCandCombPIDGenPtVsGenYXiMinus);
1086   }
1087
1088   if (!f2dHistAsMCandCombPIDGenPtVsGenYXiPlus) {
1089      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);
1090      fListHistCascade->Add(f2dHistAsMCandCombPIDGenPtVsGenYXiPlus);
1091   }
1092
1093   if (!f2dHistAsMCandCombPIDGenPtVsGenYOmegaMinus) {
1094      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);
1095      fListHistCascade->Add(f2dHistAsMCandCombPIDGenPtVsGenYOmegaMinus);
1096   }
1097
1098   if (!f2dHistAsMCandCombPIDGenPtVsGenYOmegaPlus) {
1099      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);
1100      fListHistCascade->Add(f2dHistAsMCandCombPIDGenPtVsGenYOmegaPlus);
1101   }
1102
1103
1104                 // - Generated Pt Vs Generated Y, for the cascade candidates associated with MC
1105
1106   if (!f2dHistAsMCGenPtVsGenYXiMinus) {
1107         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);
1108           fListHistCascade->Add(f2dHistAsMCGenPtVsGenYXiMinus );
1109   }
1110
1111   if (!f2dHistAsMCGenPtVsGenYXiPlus) {
1112           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);
1113           fListHistCascade->Add(f2dHistAsMCGenPtVsGenYXiPlus );
1114   }
1115
1116   if (!f2dHistAsMCGenPtVsGenYOmegaMinus) {
1117           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);
1118           fListHistCascade->Add(f2dHistAsMCGenPtVsGenYOmegaMinus );
1119   }
1120
1121   if (!f2dHistAsMCGenPtVsGenYOmegaPlus) {
1122           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);
1123           fListHistCascade->Add(f2dHistAsMCGenPtVsGenYOmegaPlus );
1124   }
1125
1126
1127                 // - Generated Eta of the the cascade candidates associated with MC
1128   if (!fHistAsMCGenEtaXiMinus) {
1129           fHistAsMCGenEtaXiMinus = new TH1F("fHistAsMCGenEtaXiMinus", "#eta of gen. #Xi^{-} (associated);#eta;Number of Casc", 100, -5, 5);
1130           fListHistCascade->Add( fHistAsMCGenEtaXiMinus );
1131   }
1132
1133   if (!fHistAsMCGenEtaXiPlus) {
1134           fHistAsMCGenEtaXiPlus = new TH1F("fHistAsMCGenEtaXiPlus", "#eta of gen. #Xi^{+} (associated);#eta;Number of Casc", 100, -5, 5);
1135           fListHistCascade->Add( fHistAsMCGenEtaXiPlus );
1136   }
1137
1138   if (!fHistAsMCGenEtaOmegaMinus) {
1139           fHistAsMCGenEtaOmegaMinus = new TH1F("fHistAsMCGenEtaOmegaMinus", "#eta of gen. #Omega^{-} (associated);#eta;Number of Casc", 100, -5, 5);
1140           fListHistCascade->Add( fHistAsMCGenEtaOmegaMinus );
1141   }
1142
1143   if (!fHistAsMCGenEtaOmegaPlus) {
1144           fHistAsMCGenEtaOmegaPlus = new TH1F("fHistAsMCGenEtaOmegaPlus", "#eta of gen. #Omega^{+} (associated);#eta;Number of Casc", 100, -5, 5);
1145           fListHistCascade->Add( fHistAsMCGenEtaOmegaPlus );
1146   }
1147
1148
1149
1150                 // - Resolution in Pt as function of generated Pt
1151
1152   if(! f2dHistAsMCResPtXiMinus) {
1153           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);
1154           fListHistCascade->Add(f2dHistAsMCResPtXiMinus);
1155   }
1156
1157   if(! f2dHistAsMCResPtXiPlus) {
1158           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);
1159           fListHistCascade->Add(f2dHistAsMCResPtXiPlus);
1160   }
1161
1162   if(! f2dHistAsMCResPtOmegaMinus) {
1163           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);
1164           fListHistCascade->Add(f2dHistAsMCResPtOmegaMinus);
1165   }
1166
1167   if(! f2dHistAsMCResPtOmegaPlus) {
1168           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);
1169           fListHistCascade->Add(f2dHistAsMCResPtOmegaPlus);
1170   }
1171
1172                 // - Resolution in R(2D) as function of generated R
1173
1174   if(! f2dHistAsMCResRXiMinus) {
1175           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);
1176           fListHistCascade->Add(f2dHistAsMCResRXiMinus);
1177   }
1178
1179   if(! f2dHistAsMCResRXiPlus) {
1180           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);
1181           fListHistCascade->Add(f2dHistAsMCResRXiPlus);
1182   }
1183
1184   if(! f2dHistAsMCResROmegaMinus) {
1185           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);
1186           fListHistCascade->Add(f2dHistAsMCResROmegaMinus);
1187   }
1188
1189   if(! f2dHistAsMCResROmegaPlus) {
1190           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);
1191           fListHistCascade->Add(f2dHistAsMCResROmegaPlus);
1192   }
1193
1194                 // - Resolution in phi as function of generated Pt
1195
1196   if(! f2dHistAsMCResPhiXiMinus) {
1197           f2dHistAsMCResPhiXiMinus = new TH2F( "f2dHistAsMCResPhiXiMinus", "Resolution in #phi for #Xi^{-}; Pt_{MC} (GeV/c); #phi(MC) - #phi(reco)   (deg)", 200, 0., 10., 60, -30., 30.);
1198           fListHistCascade->Add(f2dHistAsMCResPhiXiMinus);
1199   }
1200
1201   if(! f2dHistAsMCResPhiXiPlus) {
1202           f2dHistAsMCResPhiXiPlus = new TH2F( "f2dHistAsMCResPhiXiPlus", "Resolution in #phi for #Xi^{+}; Pt_{MC} (GeV/c); #phi(MC) - #phi(reco)   (deg)", 200, 0., 10., 60, -30., 30.);
1203           fListHistCascade->Add(f2dHistAsMCResPhiXiPlus);
1204   }
1205
1206   if(! f2dHistAsMCResPhiOmegaMinus) {
1207           f2dHistAsMCResPhiOmegaMinus = new TH2F( "f2dHistAsMCResPhiOmegaMinus", "Resolution in #phi for #Omega^{-}; Pt_{MC} (GeV/c); #phi(MC) - #phi(reco)   (deg)", 200, 0., 10., 60, -30., 30.);
1208           fListHistCascade->Add(f2dHistAsMCResPhiOmegaMinus);
1209   }
1210
1211   if(! f2dHistAsMCResPhiOmegaPlus) {
1212           f2dHistAsMCResPhiOmegaPlus = new TH2F( "f2dHistAsMCResPhiOmegaPlus", "Resolution in #phi for #Omega^{+}; Pt_{MC} (GeV/c); #phi(MC) - #phi(reco)   (deg)", 200, 0., 10., 60, -30., 30.);
1213           fListHistCascade->Add(f2dHistAsMCResPhiOmegaPlus);
1214   }
1215
1216
1217                 // - PID container
1218 if(! fCFContCascadePIDAsXiMinus&&fUseCFCont)  {
1219   const Int_t  lNbSteps      =  7 ;
1220   const Int_t  lNbVariables  =  4 ;
1221
1222   //array for the number of bins in each dimension :
1223   Int_t lNbBinsPerVar[4] = {0};
1224   lNbBinsPerVar[0] = 100;
1225   lNbBinsPerVar[1] = 800;
1226   lNbBinsPerVar[2] = 22;
1227   if(fCollidingSystems) lNbBinsPerVar[3] = 11;
1228   else lNbBinsPerVar[3] = 100;
1229
1230   fCFContCascadePIDAsXiMinus = new AliCFContainer("fCFContCascadePIDAsXiMinus","Pt_{cascade} Vs M_{#Xi^{-} candidates} Vs Y_{#Xi}", lNbSteps, lNbVariables, lNbBinsPerVar );
1231
1232   //setting the bin limits
1233   fCFContCascadePIDAsXiMinus->SetBinLimits(0,   0.0  ,  10.0 ); // Pt(Cascade)
1234   fCFContCascadePIDAsXiMinus->SetBinLimits(1,   1.2  ,   2.0 ); // Xi Effective mass
1235   fCFContCascadePIDAsXiMinus->SetBinLimits(2,  -1.1  ,   1.1 ); // Rapidity
1236   if(fCollidingSystems) {
1237     Double_t *lBinLim3  = new Double_t[ lNbBinsPerVar[3]+1 ];
1238     for(Int_t i=3; i< lNbBinsPerVar[3]+1;i++)   lBinLim3[i]  = (Double_t)(i-1)*10.;
1239     lBinLim3[0] = 0.0;
1240     lBinLim3[1] = 5.0;
1241     lBinLim3[2] = 10.0;
1242     fCFContCascadePIDAsXiMinus->SetBinLimits(3,  lBinLim3 );       // Centrality
1243   } else
1244     fCFContCascadePIDAsXiMinus->SetBinLimits(3, 0.0, 250.0  );     // SPD tracklet multiplicity
1245
1246   // Setting the step title : one per PID case
1247   fCFContCascadePIDAsXiMinus->SetStepTitle(0, "No PID");
1248   fCFContCascadePIDAsXiMinus->SetStepTitle(1, "TPC PID / 4-#sigma cut on Bachelor track");
1249   fCFContCascadePIDAsXiMinus->SetStepTitle(2, "TPC PID / 4-#sigma cut on Bachelor+Baryon tracks");
1250   fCFContCascadePIDAsXiMinus->SetStepTitle(3, "TPC PID / 4-#sigma cut on Bachelor+Baryon+Meson tracks");
1251   fCFContCascadePIDAsXiMinus->SetStepTitle(4, "Comb. PID / Bachelor");
1252   fCFContCascadePIDAsXiMinus->SetStepTitle(5, "Comb. PID / Bachelor+Baryon");
1253   fCFContCascadePIDAsXiMinus->SetStepTitle(6, "Comb. PID / Bachelor+Baryon+Meson");
1254
1255   // Setting the variable title, per axis
1256   fCFContCascadePIDAsXiMinus->SetVarTitle(0, "Pt_{cascade} (GeV/c)");
1257   fCFContCascadePIDAsXiMinus->SetVarTitle(1, "M( #Lambda , #pi^{-} ) (GeV/c^{2})");
1258   fCFContCascadePIDAsXiMinus->SetVarTitle(2, "Y_{#Xi}");
1259   if (fCollidingSystems) fCFContCascadePIDAsXiMinus->SetVarTitle(3, "Centrality");
1260   else fCFContCascadePIDAsXiMinus->SetVarTitle(3, "SPD tracklets Multiplicity");
1261
1262   fListHistCascade->Add(fCFContCascadePIDAsXiMinus);
1263
1264 }
1265
1266 if(! fCFContCascadePIDAsXiPlus&&fUseCFCont)  {
1267   const Int_t  lNbSteps      =  7 ;
1268   const Int_t  lNbVariables  =  4 ;
1269
1270   //array for the number of bins in each dimension :
1271   Int_t lNbBinsPerVar[4] = {0};
1272   lNbBinsPerVar[0] = 100;
1273   lNbBinsPerVar[1] = 800;
1274   lNbBinsPerVar[2] = 22;
1275   if(fCollidingSystems) lNbBinsPerVar[3] = 11;
1276   else lNbBinsPerVar[3] = 100;
1277
1278   fCFContCascadePIDAsXiPlus = new AliCFContainer("fCFContCascadePIDAsXiPlus","Pt_{cascade} Vs M_{#Xi^{+} candidates} Vs Y_{#Xi}", lNbSteps, lNbVariables, lNbBinsPerVar );
1279
1280
1281   //setting the bin limits (valid  for v4-18-10-AN)
1282   fCFContCascadePIDAsXiPlus->SetBinLimits(0,   0.0  ,  10.0 );  // Pt(Cascade)
1283   fCFContCascadePIDAsXiPlus->SetBinLimits(1,   1.2  ,   2.0 );  // Xi Effective mass
1284   fCFContCascadePIDAsXiPlus->SetBinLimits(2,  -1.1  ,   1.1 );  // Rapidity
1285   if(fCollidingSystems) {
1286     Double_t *lBinLim3  = new Double_t[ lNbBinsPerVar[3]+1 ];
1287     for(Int_t i=3; i< lNbBinsPerVar[3]+1;i++)   lBinLim3[i]  = (Double_t)(i-1)*10.;
1288     lBinLim3[0] = 0.0;
1289     lBinLim3[1] = 5.0;
1290     lBinLim3[2] = 10.0;
1291         fCFContCascadePIDAsXiPlus->SetBinLimits(3,lBinLim3);     // Centrality
1292   } else
1293         fCFContCascadePIDAsXiPlus->SetBinLimits(3, 0.0, 250.0  );     // SPD tracklets Multiplicity
1294
1295   // Setting the step title : one per PID case
1296   fCFContCascadePIDAsXiPlus->SetStepTitle(0, "No PID");
1297   fCFContCascadePIDAsXiPlus->SetStepTitle(1, "TPC PID / 4-#sigma cut on Bachelor track");
1298   fCFContCascadePIDAsXiPlus->SetStepTitle(2, "TPC PID / 4-#sigma cut on Bachelor+Baryon tracks");
1299   fCFContCascadePIDAsXiPlus->SetStepTitle(3, "TPC PID / 4-#sigma cut on Bachelor+Baryon+Meson tracks");
1300   fCFContCascadePIDAsXiPlus->SetStepTitle(4, "Comb. PID / Bachelor");
1301   fCFContCascadePIDAsXiPlus->SetStepTitle(5, "Comb. PID / Bachelor+Baryon");
1302   fCFContCascadePIDAsXiPlus->SetStepTitle(6, "Comb. PID / Bachelor+Baryon+Meson");
1303
1304   // Setting the variable title, per axis
1305   fCFContCascadePIDAsXiPlus->SetVarTitle(0, "Pt_{cascade} (GeV/c)");
1306   fCFContCascadePIDAsXiPlus->SetVarTitle(1, "M( #Lambda , #pi^{+} ) (GeV/c^{2})");
1307   fCFContCascadePIDAsXiPlus->SetVarTitle(2, "Y_{#Xi}");
1308   if(fCollidingSystems) fCFContCascadePIDAsXiPlus->SetVarTitle(3, "Centrality");
1309   else fCFContCascadePIDAsXiPlus->SetVarTitle(3, "SPD tracklets Multiplicity");
1310
1311   fListHistCascade->Add(fCFContCascadePIDAsXiPlus);
1312
1313 }
1314
1315
1316 if(! fCFContCascadePIDAsOmegaMinus&&fUseCFCont)  {
1317   const Int_t  lNbSteps      =  7 ;
1318   const Int_t  lNbVariables  =  4 ;
1319
1320   //array for the number of bins in each dimension :
1321   Int_t lNbBinsPerVar[4] = {0};
1322   lNbBinsPerVar[0] = 100;
1323   lNbBinsPerVar[1] = 1000;
1324   lNbBinsPerVar[2] = 22;
1325   if(fCollidingSystems) lNbBinsPerVar[3] = 11;
1326   else lNbBinsPerVar[3] = 100;
1327
1328
1329   fCFContCascadePIDAsOmegaMinus = new AliCFContainer("fCFContCascadePIDAsOmegaMinus","Pt_{cascade} Vs M_{#Omega^{-} candidates} Vs Y_{#Omega}", lNbSteps, lNbVariables, lNbBinsPerVar );
1330
1331
1332   //setting the bin limits
1333   fCFContCascadePIDAsOmegaMinus->SetBinLimits(0,   0.0  ,  10.0 );      // Pt(Cascade)
1334   fCFContCascadePIDAsOmegaMinus->SetBinLimits(1,   1.5  ,   2.5 );      // Omega Effective mass
1335   fCFContCascadePIDAsOmegaMinus->SetBinLimits(2,  -1.1  ,   1.1 );      // Rapidity
1336   if(fCollidingSystems) {
1337     Double_t *lBinLim3  = new Double_t[ lNbBinsPerVar[3]+1 ];
1338     for(Int_t i=3; i< lNbBinsPerVar[3]+1;i++)   lBinLim3[i]  = (Double_t)(i-1)*10.;
1339     lBinLim3[0] = 0.0;
1340     lBinLim3[1] = 5.0;
1341     lBinLim3[2] = 10.0;
1342         fCFContCascadePIDAsOmegaMinus->SetBinLimits(3,lBinLim3);     // Centrality
1343   } else
1344         fCFContCascadePIDAsOmegaMinus->SetBinLimits(3, 0.0, 250.0  );     // SPD tracklets multiplicity
1345
1346   // Setting the step title : one per PID case
1347   fCFContCascadePIDAsOmegaMinus->SetStepTitle(0, "No PID");
1348   fCFContCascadePIDAsOmegaMinus->SetStepTitle(1, "TPC PID / 4-#sigma cut on Bachelor track");
1349   fCFContCascadePIDAsOmegaMinus->SetStepTitle(2, "TPC PID / 4-#sigma cut on Bachelor+Baryon tracks");
1350   fCFContCascadePIDAsOmegaMinus->SetStepTitle(3, "TPC PID / 4-#sigma cut on Bachelor+Baryon+Meson tracks");
1351   fCFContCascadePIDAsOmegaMinus->SetStepTitle(4, "Comb. PID / Bachelor");
1352   fCFContCascadePIDAsOmegaMinus->SetStepTitle(5, "Comb. PID / Bachelor+Baryon");
1353   fCFContCascadePIDAsOmegaMinus->SetStepTitle(6, "Comb. PID / Bachelor+Baryon+Meson");
1354
1355   // Setting the variable title, per axis
1356   fCFContCascadePIDAsOmegaMinus->SetVarTitle(0, "Pt_{cascade} (GeV/c)");
1357   fCFContCascadePIDAsOmegaMinus->SetVarTitle(1, "M( #Lambda , K^{-} ) (GeV/c^{2})");
1358   fCFContCascadePIDAsOmegaMinus->SetVarTitle(2, "Y_{#Omega}");
1359   if(fCollidingSystems) fCFContCascadePIDAsOmegaMinus->SetVarTitle(3, "Centrality");
1360   else fCFContCascadePIDAsOmegaMinus->SetVarTitle(3, "SPD tracklet multiplicity");
1361
1362
1363   fListHistCascade->Add(fCFContCascadePIDAsOmegaMinus);
1364
1365 }
1366
1367 if(! fCFContCascadePIDAsOmegaPlus&&fUseCFCont)  {
1368   const Int_t  lNbSteps      =  7 ;
1369   const Int_t  lNbVariables  =  4 ;
1370
1371   //array for the number of bins in each dimension :
1372   Int_t lNbBinsPerVar[4]= {0};
1373   lNbBinsPerVar[0] = 100;
1374   lNbBinsPerVar[1] = 1000;
1375   lNbBinsPerVar[2] = 22;
1376   if(fCollidingSystems) lNbBinsPerVar[3] = 11;
1377   else lNbBinsPerVar[3] = 100;
1378
1379   fCFContCascadePIDAsOmegaPlus = new AliCFContainer("fCFContCascadePIDAsOmegaPlus","Pt_{cascade} Vs M_{#Omega^{+} candidates} Vs Y_{#Omega}", lNbSteps, lNbVariables, lNbBinsPerVar );
1380
1381
1382   //setting the bin limits
1383   fCFContCascadePIDAsOmegaPlus->SetBinLimits(0,   0.0  ,  10.0 );       // Pt(Cascade)
1384   fCFContCascadePIDAsOmegaPlus->SetBinLimits(1,   1.5  ,   2.5 );       // Omega Effective mass
1385   fCFContCascadePIDAsOmegaPlus->SetBinLimits(2,  -1.1  ,   1.1 );       // Rapidity
1386   if(fCollidingSystems) {
1387     Double_t *lBinLim3  = new Double_t[ lNbBinsPerVar[3]+1 ];
1388     for(Int_t i=3; i< lNbBinsPerVar[3]+1;i++)   lBinLim3[i]  = (Double_t)(i-1)*10.;
1389     lBinLim3[0] = 0.0;
1390     lBinLim3[1] = 5.0;
1391     lBinLim3[2] = 10.0;
1392         fCFContCascadePIDAsOmegaPlus->SetBinLimits(3,lBinLim3);     // Centrality
1393   } else
1394         fCFContCascadePIDAsOmegaPlus->SetBinLimits(3, 0.0, 250.0  );     // SPD tracklet multiplicity
1395
1396   // Setting the step title : one per PID case
1397   fCFContCascadePIDAsOmegaPlus->SetStepTitle(0, "No PID");
1398   fCFContCascadePIDAsOmegaPlus->SetStepTitle(1, "TPC PID / 4-#sigma cut on Bachelor track");
1399   fCFContCascadePIDAsOmegaPlus->SetStepTitle(2, "TPC PID / 4-#sigma cut on Bachelor+Baryon tracks");
1400   fCFContCascadePIDAsOmegaPlus->SetStepTitle(3, "TPC PID / 4-#sigma cut on Bachelor+Baryon+Meson tracks");
1401   fCFContCascadePIDAsOmegaPlus->SetStepTitle(4, "Comb. PID / Bachelor");
1402   fCFContCascadePIDAsOmegaPlus->SetStepTitle(5, "Comb. PID / Bachelor+Baryon");
1403   fCFContCascadePIDAsOmegaPlus->SetStepTitle(6, "Comb. PID / Bachelor+Baryon+Meson");
1404
1405   // Setting the variable title, per axis
1406   fCFContCascadePIDAsOmegaPlus->SetVarTitle(0, "Pt_{cascade} (GeV/c)");
1407   fCFContCascadePIDAsOmegaPlus->SetVarTitle(1, "M( #Lambda , K^{+} ) (GeV/c^{2})");
1408   fCFContCascadePIDAsOmegaPlus->SetVarTitle(2, "Y_{#Omega}");
1409   if(fCollidingSystems) fCFContCascadePIDAsOmegaPlus->SetVarTitle(3, "Centrality");
1410   else fCFContCascadePIDAsOmegaPlus->SetVarTitle(3,"SPD tracklet multiplicity");
1411
1412   fListHistCascade->Add(fCFContCascadePIDAsOmegaPlus);
1413
1414 }
1415
1416 // Part 3 : Towards the optimisation of topological selections -------
1417 if(! fCFContAsCascadeCuts&&fUseCFCont){
1418
1419         // Container meant to store all the relevant distributions corresponding to the cut variables.
1420         //          - FIXME optimize number of bins
1421         //          - NB overflow/underflow of variables on which we want to cut later should be 0!!!
1422
1423   const Int_t  lNbSteps      =  4 ;
1424   const Int_t  lNbVariables  =  19 ;
1425
1426   //array for the number of bins in each dimension :
1427   Int_t lNbBinsPerVar[lNbVariables] = {0};
1428   lNbBinsPerVar[0]  = 100;
1429   lNbBinsPerVar[1]  = 126;
1430   lNbBinsPerVar[2]  = 100;
1431   lNbBinsPerVar[3]  = 221;
1432   lNbBinsPerVar[4]  = 30;
1433   lNbBinsPerVar[5]  = 50;
1434
1435   lNbBinsPerVar[6]  = 100;
1436   lNbBinsPerVar[7]  = 43;
1437   lNbBinsPerVar[8]  = 101;
1438   lNbBinsPerVar[9]  = 26;
1439   lNbBinsPerVar[10] = 26;
1440
1441   lNbBinsPerVar[11] = 150; // 2-MeV/c2 bins
1442   lNbBinsPerVar[12] = 120; // 2-MeV/c2 bins
1443
1444   lNbBinsPerVar[13] = 100;
1445
1446   lNbBinsPerVar[14] = 44; // 0.05 in rapidity units
1447   lNbBinsPerVar[15] = 44; // 0.05 in rapidity units
1448
1449   lNbBinsPerVar[16] = 20;
1450
1451
1452   if(fCollidingSystems) lNbBinsPerVar[17] = 11;
1453   else lNbBinsPerVar[17] = 100;
1454   lNbBinsPerVar[18] = 100;
1455
1456
1457   fCFContAsCascadeCuts = new AliCFContainer("fCFContAsCascadeCuts","Cut Container for Asso. Cascades", lNbSteps, lNbVariables, lNbBinsPerVar );
1458
1459   //0
1460   fCFContAsCascadeCuts->SetBinLimits(0, 0., 2.);                 // DcaXiDaughters : 0.0 to 2.0
1461   //1
1462    Double_t *lBinLim1  = new Double_t[ lNbBinsPerVar[1]+1 ];
1463    for(Int_t i=0; i< lNbBinsPerVar[1];i++)   lBinLim1[i]  = (Double_t)0.0   + (5.  - 0.0 )/(lNbBinsPerVar[1]-1)  * (Double_t)i ;
1464    lBinLim1[ lNbBinsPerVar[1]  ] = 100.0;
1465    fCFContAsCascadeCuts -> SetBinLimits(1,  lBinLim1 );
1466   delete [] lBinLim1;                                            // DcaBachToPrimVertexXi : 0.0 to 0.5
1467   //2
1468   fCFContAsCascadeCuts->SetBinLimits(2, .99, 1.);                // XiCosineOfPointingAngle : 0.99 to 1.0
1469   //3
1470   Double_t *lBinLim3  = new Double_t[ lNbBinsPerVar[3]+1 ];
1471   for(Int_t i=0; i< lNbBinsPerVar[3];i++)   lBinLim3[i] = (Double_t)0.0 + (5.  - 0.0 )/(lNbBinsPerVar[3]-1) * (Double_t)i ;
1472   lBinLim3[ lNbBinsPerVar[3]  ] = 110.0;
1473   fCFContAsCascadeCuts -> SetBinLimits(3,  lBinLim3 );            // XiRadius : 0.0 to 4.0
1474   delete [] lBinLim3;
1475
1476   //4
1477   fCFContAsCascadeCuts->SetBinLimits(4, 1.1, 1.13);               // InvMassLambdaAsCascDghter
1478   //5
1479   fCFContAsCascadeCuts->SetBinLimits(5, 0., 2.);                  // DcaV0DaughtersXi : 0.0 to 2.0
1480   //6
1481   fCFContAsCascadeCuts->SetBinLimits(6, 0.99, 1.);                // V0CosineOfPointingAngleXi : 0.98 to 1.0
1482   //7
1483   Double_t *lBinLim7  = new Double_t[ lNbBinsPerVar[7]+1 ];
1484   for(Int_t i=0; i< lNbBinsPerVar[7]-2;i++)   lBinLim7[i]  = (Double_t)0.0   + (20.  - 0.0 )/(lNbBinsPerVar[7]-3)  * (Double_t)i ;
1485   lBinLim7[ lNbBinsPerVar[7]-2] = 100.0;
1486   lBinLim7[ lNbBinsPerVar[7]-1] = 200.0;
1487   lBinLim7[ lNbBinsPerVar[7]] = 1000.0;
1488   fCFContAsCascadeCuts -> SetBinLimits(7,  lBinLim7 );
1489   delete [] lBinLim7;                                             // V0RadiusXi : 0.0 to 20.0
1490   //8
1491   Double_t *lBinLim8  = new Double_t[ lNbBinsPerVar[8]+1 ];
1492         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 ;
1493         lBinLim8[ lNbBinsPerVar[8]  ] = 100.0;
1494   fCFContAsCascadeCuts -> SetBinLimits(8,  lBinLim8 );
1495   delete [] lBinLim8;                                             // DcaV0ToPrimVertexXi : 0. to 0.4
1496   //9
1497   Double_t *lBinLim9  = new Double_t[ lNbBinsPerVar[9]+1 ];
1498   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 ;
1499   lBinLim9[ lNbBinsPerVar[9]  ] = 100.0;
1500   fCFContAsCascadeCuts -> SetBinLimits(9,  lBinLim9 );
1501   delete [] lBinLim9;                                             // DcaPosToPrimVertexXi : 0.0 to 0.25
1502   //10
1503   Double_t *lBinLim10  = new Double_t[ lNbBinsPerVar[10]+1 ];
1504   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 ;
1505   lBinLim10[ lNbBinsPerVar[10]  ] = 100.0;
1506   fCFContAsCascadeCuts -> SetBinLimits(10,  lBinLim10 );
1507   delete [] lBinLim10;                                            // DcaPosToPrimVertexXi : 0.0 to 0.25
1508   //11
1509   fCFContAsCascadeCuts->SetBinLimits(11, 1.25, 1.40);             // InvMassXi
1510   fCFContAsCascadeCuts->SetBinLimits(12, 1.62, 1.74);             // InvMassOmega
1511   fCFContAsCascadeCuts->SetBinLimits(13, 0.0, 10.0);              // XiTransvMom
1512   fCFContAsCascadeCuts->SetBinLimits(14, -1.1, 1.1);              // Y(Xi)
1513   fCFContAsCascadeCuts->SetBinLimits(15, -1.1, 1.1);              // Y(Omega)
1514   fCFContAsCascadeCuts->SetBinLimits(16, -10.0, 10.0);            // BestPrimaryVtxPosZ
1515   if (fCollidingSystems) {
1516     Double_t *lBinLim17  = new Double_t[ lNbBinsPerVar[17]+1 ];
1517     for(Int_t i=3; i< lNbBinsPerVar[17]+1;i++)   lBinLim17[i]  = (Double_t)(i-1)*10.;
1518     lBinLim17[0] = 0.0;
1519     lBinLim17[1] = 5.0;
1520     lBinLim17[2] = 10.0;
1521     fCFContAsCascadeCuts -> SetBinLimits(17,  lBinLim17 );       // Centrality
1522     delete [] lBinLim17;
1523     fCFContAsCascadeCuts->SetBinLimits(18, 0.0, 6000.0);         // ESD track multiplicity
1524   } else {
1525     fCFContAsCascadeCuts->SetBinLimits(17, 0.0, 250.0);          // SPDTrackletsMultiplicity
1526     fCFContAsCascadeCuts->SetBinLimits(18, 0.0, 200.0);          // ESD track multiplicity
1527   }
1528
1529   // Setting the number of steps : one for each cascade species (Xi-, Xi+ and Omega-, Omega+)
1530   fCFContAsCascadeCuts->SetStepTitle(0, "#Xi^{-} candidates associated to MC");
1531   fCFContAsCascadeCuts->SetStepTitle(1, "#bar{#Xi}^{+} candidates associated to MC");
1532   fCFContAsCascadeCuts->SetStepTitle(2, "#Omega^{-} candidates associated to MC");
1533   fCFContAsCascadeCuts->SetStepTitle(3, "#bar{#Omega}^{+} candidates associated to MC");
1534
1535   // Setting the variable title, per axis
1536   fCFContAsCascadeCuts->SetVarTitle(0,  "DCA(XiDaughters) (cm)");
1537   fCFContAsCascadeCuts->SetVarTitle(1,  "DCA(Bach/PrimVertex) (cm)");
1538   fCFContAsCascadeCuts->SetVarTitle(2,  "cos(Xi pointing angle)");
1539   fCFContAsCascadeCuts->SetVarTitle(3,  "R_{2d}(Xi decay) (cm)");
1540   fCFContAsCascadeCuts->SetVarTitle(4,  "M_{#Lambda}(As Casc Dghter) (GeV/c^{2})");
1541   fCFContAsCascadeCuts->SetVarTitle(5,  "DCA(V0 Daughters Xi) (cm)");
1542
1543   fCFContAsCascadeCuts->SetVarTitle(6,  "cos(V0 pointing Angle) in Casc");
1544   fCFContAsCascadeCuts->SetVarTitle(7,  "R_{2d}(V0 decay) (cm)");
1545   fCFContAsCascadeCuts->SetVarTitle(8,  "DCA(V0/PrimVertex) (cm)");
1546   fCFContAsCascadeCuts->SetVarTitle(9,  "DCA(Pos/PrimVertex) (cm)");
1547   fCFContAsCascadeCuts->SetVarTitle(10, "DCA(Neg/PrimVertex) (cm)");
1548
1549   fCFContAsCascadeCuts->SetVarTitle(11, "Inv. Mass(Xi) (GeV/c^{2})");
1550   fCFContAsCascadeCuts->SetVarTitle(12, "Inv. Mass(Omega) (GeV/c^{2})");
1551
1552   fCFContAsCascadeCuts->SetVarTitle(13, "Pt_{MC}(Casc.) (GeV/c)");
1553
1554   fCFContAsCascadeCuts->SetVarTitle(14, "Y_{MC}(Xi)");
1555   fCFContAsCascadeCuts->SetVarTitle(15, "Y_{MC}(Omega)");
1556
1557   fCFContAsCascadeCuts->SetVarTitle(16, "Z-position(BestPrimVtx) (cm)");
1558
1559   if (fCollidingSystems) fCFContAsCascadeCuts->SetVarTitle(17, "Centrality");
1560   else fCFContAsCascadeCuts->SetVarTitle(17, "SPD tracklets Multiplicity");
1561   fCFContAsCascadeCuts->SetVarTitle(18, "ESD track multiplicity");
1562
1563   fListHistCascade->Add(fCFContAsCascadeCuts);
1564 }
1565
1566   fV0Ampl = new TH1F("fV0Ampl","",500,0.,30000);
1567   fListHistCascade->Add(fV0Ampl);
1568
1569
1570 PostData(1, fListHistCascade);
1571 }// end CreateOutputObjects
1572
1573
1574 //________________________________________________________________________
1575 void AliAnalysisTaskCheckPerformanceCascadePbPb::UserExec(Option_t *)
1576 {
1577
1578   // Main loop
1579   // Called for each event
1580
1581         AliESDEvent *lESDevent = 0x0;
1582         AliAODEvent *lAODevent = 0x0;
1583         AliMCEvent  *lMCevent  = 0x0;
1584         AliStack    *lMCstack  = 0x0;
1585         Int_t ncascades = -1;
1586
1587
1588   // Connect to the InputEvent
1589   // After these lines, we should have an ESD/AOD event + the number of cascades in it.
1590
1591         if (fAnalysisType == "ESD") {
1592                 lESDevent = dynamic_cast<AliESDEvent*>( InputEvent() );
1593                 if (!lESDevent) {
1594                         Printf("ERROR: lESDevent not available \n");
1595                         cout << "Name of the file with pb :" <<  CurrentFileName() << endl;  // or AliAnalysisTaskSE::CurrentFileName()
1596                         return;
1597                 }
1598         } else if (fAnalysisType == "AOD") {
1599                 lAODevent = dynamic_cast<AliAODEvent*>( InputEvent() );
1600                 if (!lAODevent) {
1601                         Printf("ERROR: lAODevent not available \n");
1602                         cout << "Name of the file with pb :" <<  CurrentFileName() << endl;
1603                         return;
1604                 }
1605         }
1606
1607
1608         lMCevent = MCEvent();
1609         if (!lMCevent) {
1610                 Printf("ERROR: Could not retrieve MC event \n");
1611                 cout << "Name of the file with pb :" <<  CurrentFileName() << endl;
1612                 return;
1613         }
1614
1615         lMCstack = lMCevent->Stack();
1616         if (!lMCstack) {
1617                 Printf("ERROR: Could not retrieve MC stack \n");
1618                 cout << "Name of the file with pb :" <<  CurrentFileName() << endl;
1619                 return;
1620
1621         }
1622
1623
1624    // Temporary way : AOD awareness of the code to be developed  FIXME
1625    if (fAnalysisType == "AOD") return;
1626
1627   //-------------------------------------------------
1628   // 1 - Cascade vertexer (ESD)
1629         if (fkRerunV0CascVertexers) { // relaunch V0 and Cascade vertexers, not test in PbPb
1630                 if(fAnalysisType == "ESD" ){
1631 //                         lESDevent->ResetCascades();
1632 //                         lESDevent->ResetV0s();
1633 //
1634 //                         AliV0vertexer lV0vtxer;
1635 //                         AliCascadeVertexer lCascVtxer;
1636 //
1637 //                         lV0vtxer.SetDefaultCuts(fV0Sels);
1638 //                         lCascVtxer.SetDefaultCuts(fCascSels);
1639 //
1640 //                         lV0vtxer.Tracks2V0vertices(lESDevent);
1641 //                         lCascVtxer.V0sTracks2CascadeVertices(lESDevent);
1642                 }
1643         }
1644
1645
1646   //------------------------------------------------
1647   // 2 - Preparing the general info about of the event = prim. Vtx + magnetic field (ESD)
1648
1649
1650 //      if(fAnalysisType == "ESD" ){
1651
1652         // Magnetic field
1653                 const Double_t lMagneticField = lESDevent->GetMagneticField( );
1654
1655         // Prim vertex
1656                 const AliESDVertex *lPrimaryTrackingVtx = lESDevent->GetPrimaryVertexTracks();  // get the vtx stored in ESD found with tracks
1657                 const AliESDVertex *lPrimarySPDVtx      = lESDevent->GetPrimaryVertexSPD();     // get the vtx stored in ESD found with SPD tracklets
1658
1659                 const AliESDVertex *lPrimaryBestVtx     = lESDevent->GetPrimaryVertex();
1660                         // get the best primary vertex available for the event
1661                         // As done in AliCascadeVertexer, we keep the one which is the best one available.
1662                         // between : Tracking vertex > SPD vertex > TPC vertex > default SPD vertex
1663                         Double_t lBestPrimaryVtxPos[3]   = {-100.0, -100.0, -100.0};
1664                 lPrimaryBestVtx->GetXYZ( lBestPrimaryVtxPos );
1665
1666                 // FIXME : quality cut on the z-position of the prim vertex.
1667                 if(fkQualityCutZprimVtxPos) {
1668                         if(TMath::Abs(lBestPrimaryVtxPos[2]) > fVtxRange ) {
1669                                 AliWarning("Pb / | Z position of Best Prim Vtx | > 10.0 cm ... return !");
1670                                 PostData(1, fListHistCascade);
1671                                 return;
1672                         }
1673                 }
1674                 // FIXME : quality selection regarding pile-up rejection
1675                 if(fkRejectEventPileUp) {
1676                         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
1677                                 AliWarning("Pb / Event tagged as pile-up by SPD... return !");
1678                                 PostData(1, fListHistCascade);
1679                                 return;
1680                         }
1681                 }
1682                 // FIXME : remove TPC-only primary vertex : retain only events with tracking + SPD vertex
1683                 if(fkQualityCutNoTPConlyPrimVtx) {
1684                         if (!lPrimarySPDVtx->GetStatus() && !lPrimaryTrackingVtx->GetStatus() ){
1685                                 AliWarning("Pb / No SPD prim. vertex nor prim. Tracking vertex ... return !");
1686                                 PostData(1, fListHistCascade);
1687                                 return;
1688                         }
1689                 }
1690
1691   // Centrality determination
1692   AliESDVZERO* esdV0 = lESDevent->GetVZEROData();
1693   Float_t multV0A=esdV0->GetMTotV0A();
1694   Float_t multV0C=esdV0->GetMTotV0C();
1695
1696   AliCentrality *centrality = lESDevent->GetCentrality();
1697   //  Printf("Centrality percentile V0M for this event %f)\n",  centrality->GetCentralityPercentile("V0M"));
1698   Float_t lcentrality = centrality->GetCentralityPercentile(fCentrEstimator.Data());
1699   if (lcentrality==100.) lcentrality=99.9;
1700 /*  if (lcentrality==-1||lcentrality>=90.) {
1701     PostData(1, fListHistCascade);
1702     return;
1703   }
1704 */
1705 /*  if (!(centrality->IsEventInCentralityClass(fCentrLowLim,fCentrUpLim,fCentrEstimator.Data()))) {
1706     PostData(1, fListHistCascade);
1707     return;
1708   }
1709 */
1710   fV0Ampl->Fill(multV0A+multV0C);
1711
1712 //      }// if ESD
1713
1714
1715   //    cout << "Name of the accessed file :" <<  fInputHandler->GetTree()->GetCurrentFile()->GetName() << endl;
1716
1717   //    cout << "Tree characteristics ..." << endl;
1718   //    fInputHandler->GetTree()->Print("toponly");
1719   //    fInputHandler->GetTree()->GetBranch("PrimaryVertex")->Print();
1720   //    fInputHandler->GetTree()->GetBranch("SPDVertex")->Print();
1721
1722
1723
1724   // ---------------------------------------------------------------
1725   // - Initialisation of the part dedicated to cascade vertices
1726
1727   if(fAnalysisType == "ESD")            ncascades = lESDevent->GetNumberOfCascades();
1728   else if(fAnalysisType == "AOD")       ncascades = lAODevent->GetNumberOfCascades();
1729
1730
1731   Int_t   nNumberOfMCPrimaries       = -1;
1732   Int_t   nMCPrimariesInAcceptance   =  0;
1733   Int_t   nTrackPrimaryMultiplicity  = -1;
1734
1735
1736 //        nNumberOfMCPrimaries       = lMCstack->GetNprimary();
1737         nNumberOfMCPrimaries = lMCstack->GetNtrack(); // MN: this stack->GetNtrack(); has to be used because in HIJING decay products of D and B mesons are also primaries and produced in HIJING during transport
1738
1739         nTrackPrimaryMultiplicity  = fESDtrackCuts->CountAcceptedTracks(lESDevent);
1740
1741         fHistEvtsInCentralityBinsvsNtracks->Fill(lcentrality,nTrackPrimaryMultiplicity);
1742
1743         if(nNumberOfMCPrimaries < 1) return;
1744
1745         fHistMCTrackMultiplicity->Fill( nNumberOfMCPrimaries );  // MN: neutral particles included and also not physical ones
1746
1747 //_____________________________________________________________________________
1748 // Part 1 - Loop over the MC primaries
1749
1750     for (Int_t iCurrentLabelStack = 0; iCurrentLabelStack < nNumberOfMCPrimaries; iCurrentLabelStack++) {// This is the beginning of the loop on primaries
1751
1752         TParticle* lCurrentParticle = 0x0;
1753                    lCurrentParticle = lMCstack->Particle( iCurrentLabelStack );
1754         if(!lCurrentParticle) {
1755                 Printf("MC Primary loop %d - MC TParticle pointer to current stack particle = 0x0 ! Skip ...\n", iCurrentLabelStack );
1756                 continue;
1757         }
1758
1759         if (!lMCstack->IsPhysicalPrimary(iCurrentLabelStack)) continue;
1760
1761
1762         if( TMath::Abs( lCurrentParticle->Eta() ) > 0.8 ) continue;
1763         nMCPrimariesInAcceptance++;
1764     }
1765
1766     f2dHistRecoMultVsMCMult->Fill( nTrackPrimaryMultiplicity, nMCPrimariesInAcceptance );  // MN: neutral are included
1767
1768
1769    // For proton
1770
1771 /*   for (Int_t iCurrentLabelStack = 0; iCurrentLabelStack < nNumberOfMCPrimaries; iCurrentLabelStack++)
1772         {// This is the begining of the loop on primaries, for protons
1773
1774         TParticle* lCurrentParticle = lMCstack->Particle( iCurrentLabelStack );
1775         if(!lCurrentParticle){
1776                 Printf("Proton loop %d - MC TParticle pointer to current stack particle = 0x0 ! Skip ...\n", iCurrentLabelStack );
1777                 continue;
1778
1779         }
1780
1781         if( lCurrentParticle->GetPdgCode() == 2212 )
1782                 fHistEtaGenProton->Fill( lCurrentParticle->Eta() );
1783
1784         if( lCurrentParticle->GetPdgCode() == -2212 )
1785                 fHistEtaGenAntiProton->Fill( lCurrentParticle->Eta() );
1786         }// end loop over primary proton
1787 */
1788
1789
1790
1791 //_____________________________________________________________________________
1792 // Part 2 - Loop over the different types of GENERATED cascades (Xi-+, Omega-+)
1793
1794         // - Initialisation of useful local variables
1795
1796         Int_t lPdgCodeCasc            = 0;
1797         Int_t lPdgCodeBach            = 0;
1798         Int_t lPdgCodeLambda          = 0;
1799         Int_t lPdgCodeDghtMesV0       = 0;
1800         Int_t lPdgCodeDghtBarV0       = 0;
1801
1802
1803         TH1F *lHistEtaGenCasc         = 0;
1804         TH3D *l3dHistGenPtVsGenYGenvsCent = 0;
1805         TH3D *l3dHistGenPtVsGenYGenvsNtracks = 0;
1806         TH1F *lHistThetaGenCasc       = 0;
1807         TH2D *l2dHistGenPtVsGenYFdbl  = 0;
1808         TH1F *lHistThetaLambda        = 0;
1809         TH1F *lHistThetaBach          = 0;
1810         TH1F *lHistThetaBarDghter     = 0;
1811         TH1F *lHistThetaMesDghter     = 0;
1812         TH1F *lHistPtBach             = 0;
1813         TH1F *lHistPtBarDghter        = 0;
1814         TH1F *lHistPtMesDghter        = 0;
1815         Int_t ncascperev = 0;
1816         Int_t ncascperevtot = 0;
1817
1818
1819 for (Int_t iCascType = 1; iCascType < 5; iCascType++) {
1820   ncascperev = 0;
1821   ncascperevtot = 0;
1822
1823   switch (iCascType) {
1824     case 1: // Xi-
1825          lPdgCodeCasc       =   3312;  //Xi-
1826          lPdgCodeBach       =   -211;  //Pi-
1827          lPdgCodeLambda     =   3122;  //Lambda0
1828          lPdgCodeDghtMesV0  =   -211;  //Pi-
1829          lPdgCodeDghtBarV0  =   2212;  //Proton
1830
1831                 // any Xi-
1832          lHistEtaGenCasc        = fHistEtaGenCascXiMinus;
1833          l3dHistGenPtVsGenYGenvsCent = f3dHistGenPtVsGenYGenvsCentXiMinus;
1834          l3dHistGenPtVsGenYGenvsNtracks = f3dHistGenPtVsGenYGenvsNtracksXiMinus;
1835                 // cascades generated within acceptance (cut in pt + theta)
1836          lHistThetaGenCasc      = fHistThetaGenCascXiMinus;
1837          l2dHistGenPtVsGenYFdbl = f2dHistGenPtVsGenYFdblXiMinus;
1838          lHistThetaLambda       = fHistThetaLambdaXiMinus;
1839          lHistThetaBach         = fHistThetaBachXiMinus;
1840          lHistThetaBarDghter    = fHistThetaBarDghterXiMinus;
1841          lHistThetaMesDghter    = fHistThetaMesDghterXiMinus;
1842          lHistPtBach            = fHistPtBachXiMinus;
1843          lHistPtBarDghter       = fHistPtBarDghterXiMinus;
1844          lHistPtMesDghter       = fHistPtMesDghterXiMinus;
1845         break;
1846
1847     case 2: // Xi+
1848          lPdgCodeCasc        =  -3312;  //Xi+
1849          lPdgCodeBach        =    211;  //Pi+
1850          lPdgCodeLambda      =  -3122;  //AntiLambda0
1851          lPdgCodeDghtMesV0   =    211;  //Pi+
1852          lPdgCodeDghtBarV0   =  -2212;  //AntiProton
1853
1854                 // any Xi+
1855          lHistEtaGenCasc        = fHistEtaGenCascXiPlus;
1856          //l2dHistGenPtVsGenYGen  = f2dHistGenPtVsGenYGenXiPlus;
1857          l3dHistGenPtVsGenYGenvsCent = f3dHistGenPtVsGenYGenvsCentXiPlus;
1858          l3dHistGenPtVsGenYGenvsNtracks = f3dHistGenPtVsGenYGenvsNtracksXiPlus;
1859
1860
1861                 // cascades generated within acceptance (cut in pt + theta)
1862          lHistThetaGenCasc      = fHistThetaGenCascXiPlus;
1863          l2dHistGenPtVsGenYFdbl = f2dHistGenPtVsGenYFdblXiPlus;
1864          lHistThetaLambda       = fHistThetaLambdaXiPlus;
1865          lHistThetaBach         = fHistThetaBachXiPlus;
1866          lHistThetaBarDghter    = fHistThetaBarDghterXiPlus;
1867          lHistThetaMesDghter    = fHistThetaMesDghterXiPlus;
1868          lHistPtBach            = fHistPtBachXiPlus;
1869          lHistPtBarDghter       = fHistPtBarDghterXiPlus;
1870          lHistPtMesDghter       = fHistPtMesDghterXiPlus;
1871         break;
1872
1873     case 3: // Omega-
1874          lPdgCodeCasc       =   3334;  //Omega-
1875          lPdgCodeBach       =   -321;  //K-
1876          lPdgCodeLambda     =   3122;  //Lambda0
1877          lPdgCodeDghtMesV0  =   -211;  //Pi-
1878          lPdgCodeDghtBarV0  =   2212;  //Proton
1879
1880                 // any Omega-
1881          lHistEtaGenCasc        = fHistEtaGenCascOmegaMinus;
1882 //       l2dHistGenPtVsGenYGen  = f2dHistGenPtVsGenYGenOmegaMinus;
1883          l3dHistGenPtVsGenYGenvsCent = f3dHistGenPtVsGenYGenvsCentOmegaMinus;
1884          l3dHistGenPtVsGenYGenvsNtracks = f3dHistGenPtVsGenYGenvsNtracksOmegaMinus;
1885
1886
1887                 // cascades generated within acceptance (cut in pt + theta)
1888          lHistThetaGenCasc      = fHistThetaGenCascOmegaMinus;
1889          l2dHistGenPtVsGenYFdbl = f2dHistGenPtVsGenYFdblOmegaMinus;
1890          lHistThetaLambda       = fHistThetaLambdaOmegaMinus;
1891          lHistThetaBach         = fHistThetaBachOmegaMinus;
1892          lHistThetaBarDghter    = fHistThetaBarDghterOmegaMinus;
1893          lHistThetaMesDghter    = fHistThetaMesDghterOmegaMinus;
1894          lHistPtBach            = fHistPtBachOmegaMinus;
1895          lHistPtBarDghter       = fHistPtBarDghterOmegaMinus;
1896          lHistPtMesDghter       = fHistPtMesDghterOmegaMinus;
1897         break;
1898
1899     case 4:  // Omega+
1900          lPdgCodeCasc       =  -3334;  //Omega+
1901          lPdgCodeBach       =    321;  //K+
1902          lPdgCodeLambda     =  -3122;  //AntiLambda0
1903          lPdgCodeDghtMesV0  =    211;  //Pi+
1904          lPdgCodeDghtBarV0  =  -2212;  //AntiProton
1905
1906                 // any Omega+
1907          lHistEtaGenCasc        = fHistEtaGenCascOmegaPlus;
1908          l3dHistGenPtVsGenYGenvsCent = f3dHistGenPtVsGenYGenvsCentOmegaPlus;
1909          l3dHistGenPtVsGenYGenvsNtracks = f3dHistGenPtVsGenYGenvsNtracksOmegaPlus;
1910
1911
1912                 // cascades generated within acceptance (cut in pt + theta)
1913          lHistThetaGenCasc      = fHistThetaGenCascOmegaPlus;
1914          l2dHistGenPtVsGenYFdbl = f2dHistGenPtVsGenYFdblOmegaPlus;
1915          lHistThetaLambda       = fHistThetaLambdaOmegaPlus;
1916          lHistThetaBach         = fHistThetaBachOmegaPlus;
1917          lHistThetaBarDghter    = fHistThetaBarDghterOmegaPlus;
1918          lHistThetaMesDghter    = fHistThetaMesDghterOmegaPlus;
1919          lHistPtBach            = fHistPtBachOmegaPlus;
1920          lHistPtBarDghter       = fHistPtBarDghterOmegaPlus;
1921          lHistPtMesDghter       = fHistPtMesDghterOmegaPlus;
1922         break;
1923
1924   }// end switch cascade
1925
1926
1927   for (Int_t iCurrentLabelStack = 0; iCurrentLabelStack < nNumberOfMCPrimaries; iCurrentLabelStack++) {// This is the beginning of the loop on primaries
1928
1929         TParticle* lCurrentParticle = 0x0;
1930                    lCurrentParticle = lMCstack->Particle( iCurrentLabelStack );
1931         if (!lCurrentParticle) {
1932           Printf("Cascade loop %d - MC TParticle pointer to current stack particle = 0x0 ! Skip ...\n", iCurrentLabelStack );
1933           continue;
1934         }
1935         if (!lMCstack->IsPhysicalPrimary(iCurrentLabelStack)) continue;
1936         if (lCurrentParticle->GetPdgCode() == lPdgCodeCasc) {  // Here !
1937                 //cout << "Xi- within loop " << iCurrentLabelStack << "/ " << nNumberOfMCPrimaries << endl;
1938
1939                 // -  Xi level ... _____________________________________________________________
1940                 TParticle* xiMC = 0x0;
1941                            xiMC = lCurrentParticle;
1942                 if(!xiMC){
1943                         Printf("MC TParticle pointer to Cascade = 0x0 ! Skip ...");
1944                         continue;
1945                 } // redundant?
1946
1947                 // To select injected/bio particles Bool_t AliMCEvent::IsFromBGEvent(Int_t index)
1948
1949                 ncascperevtot++;
1950                 // Fill the first histos : = any generated Xi, not necessarily within the acceptance
1951                 Double_t lRapXiMC = 0.5*TMath::Log((xiMC->Energy() + xiMC->Pz()) / (xiMC->Energy() - xiMC->Pz() +1.e-13));
1952
1953                 lHistEtaGenCasc         ->Fill( xiMC->Eta() );
1954                 l3dHistGenPtVsGenYGenvsCent   ->Fill( xiMC->Pt(), lRapXiMC, lcentrality );
1955                 l3dHistGenPtVsGenYGenvsNtracks->Fill( xiMC->Pt(), lRapXiMC, nTrackPrimaryMultiplicity );
1956
1957
1958
1959                 // Check the emission of particle stays within the acceptance of the detector (cut in theta)
1960                 if (fApplyAccCut) {if( xiMC->Theta() < TMath::Pi()/4.0  ||    xiMC->Theta() > 3.0*TMath::Pi()/4.0 ) continue;}
1961                 if( xiMC->GetNDaughters() != 2) continue;
1962                 if( xiMC->GetDaughter(0) < 0 )  continue;
1963                 if( xiMC->GetDaughter(1) < 0 )  continue;
1964
1965                         TParticle* lDght0ofXi = lMCstack->Particle(  xiMC->GetDaughter(0) );
1966                         TParticle* lDght1ofXi = lMCstack->Particle(  xiMC->GetDaughter(1) );
1967
1968                 TParticle* lLambda = 0;
1969                 TParticle* lBach   = 0;
1970
1971                 // Xi - Case 1
1972                         if(     lDght0ofXi->GetPdgCode() == lPdgCodeLambda   &&  // Here !
1973                                 lDght1ofXi->GetPdgCode() == lPdgCodeBach ){      // Here !
1974
1975                                 lLambda = lDght0ofXi;
1976                                 lBach   = lDght1ofXi;
1977                         }// end if dghter 0 = Lambda and    dghter 1 = Pi-
1978
1979                 // Xi - Case 2
1980                         else if( lDght0ofXi->GetPdgCode() == lPdgCodeBach  &&      // Here !
1981                                  lDght1ofXi->GetPdgCode() == lPdgCodeLambda ){     // Here !
1982
1983                                 lBach   = lDght0ofXi;
1984                                 lLambda = lDght1ofXi;
1985                         }//  end if dghter 0 = Pi-  and   dghter 1 = Lambda
1986
1987                 // V0 otherwise - Case 3
1988                         else continue;
1989
1990                         // Check the emission of particle stays within the acceptance of the detector (cut in pt + theta)
1991                         if (fApplyAccCut) {
1992                           if( lLambda->Theta() < TMath::Pi()/4.0  ||    lLambda->Theta() > 3.0*TMath::Pi()/4.0 ) continue;
1993                           if( lBach->Theta() < TMath::Pi()/4.0    ||    lBach->Theta() > 3.0*TMath::Pi()/4.0 )   continue;
1994                           if( lBach->Pt() < 0.150 ) continue; //FIXME : maybe tuned for Xi but not for K- from Omega ...
1995                         }
1996
1997                 // -  V0 level ... _____________________________________________________________
1998                 TParticle* lDghtBarV0 = 0;
1999                 TParticle* lDghtMesV0 = 0;
2000
2001                 if( lLambda->GetNDaughters() != 2 )  continue;
2002                 if( lLambda->GetDaughter(0) < 0 )    continue;
2003                 if( lLambda->GetDaughter(1) < 0 )    continue;
2004
2005
2006                 TParticle* lDght0ofLambda = lMCstack->Particle(  lLambda->GetDaughter(0) );
2007                 TParticle* lDght1ofLambda = lMCstack->Particle(  lLambda->GetDaughter(1) );
2008
2009                 // V0 - Case 1
2010                         if(     lDght0ofLambda->GetPdgCode() == lPdgCodeDghtBarV0 &&    // Here !
2011                                 lDght1ofLambda->GetPdgCode() == lPdgCodeDghtMesV0 ){    // Here !
2012
2013                                 lDghtBarV0 = lDght0ofLambda;
2014                                 lDghtMesV0 = lDght1ofLambda;
2015                         }// end if dghter 0 = Proton  and   dghter 1 = Pi-
2016
2017                 // V0 - Case 2
2018                         else if( lDght0ofLambda->GetPdgCode() == lPdgCodeDghtMesV0  &&     // Here !
2019                                  lDght1ofLambda->GetPdgCode() == lPdgCodeDghtBarV0 ){      // Here !
2020
2021                                 lDghtMesV0 = lDght0ofLambda;
2022                                 lDghtBarV0 = lDght1ofLambda;
2023                         }//  end if dghter 0 = Pi-  and   dghter 1 = proton
2024
2025                 // V0 otherwise - Case 3
2026                         else continue;
2027
2028
2029                         // Check the emission of particle stays within the acceptance of the detector
2030                         if (fApplyAccCut) {
2031                           if( lDghtBarV0->Theta() < TMath::Pi()/4.0  ||  lDghtBarV0->Theta() > 3.0*TMath::Pi()/4.0 ) continue;
2032                           if( lDghtMesV0->Theta() < TMath::Pi()/4.0  ||  lDghtMesV0->Theta() > 3.0*TMath::Pi()/4.0 ) continue;
2033
2034                           if( lDghtBarV0->Pt() < 0.250 ) continue;
2035                           if( lDghtMesV0->Pt() < 0.150 ) continue;
2036                         }
2037
2038
2039                 // - Just to know which file is currently open : locate the file containing Xi
2040                 //cout << "Name of the file containing generated Xi :" <<  fInputHandler->GetTree()->GetCurrentFile()->GetName()
2041                 //                                                   <<  endl;
2042
2043                 Double_t lRadToDeg = 180.0/TMath::Pi();
2044
2045                 // - Filling histos ... _________________________________________________________________
2046                         lHistThetaGenCasc       ->Fill( lRadToDeg * xiMC->Theta()  );
2047                         l2dHistGenPtVsGenYFdbl  ->Fill( xiMC->Pt(), lRapXiMC );
2048
2049                         // - Fill theta histos for Lambda and Bach
2050                         lHistThetaLambda        ->Fill( lRadToDeg * lLambda->Theta() );
2051                         lHistThetaBach          ->Fill( lRadToDeg *   lBach->Theta() );
2052
2053                         // - Fill theta histos for V0 daughters
2054                         lHistThetaBarDghter     ->Fill( lRadToDeg * lDghtBarV0->Theta() );
2055                         lHistThetaMesDghter     ->Fill( lRadToDeg * lDghtMesV0->Theta() );
2056
2057                         // - Fill pt histos.
2058                         lHistPtBach             ->Fill(      lBach->Pt() );
2059                         lHistPtBarDghter        ->Fill( lDghtBarV0->Pt() );
2060                         lHistPtMesDghter        ->Fill( lDghtMesV0->Pt() );
2061                       //  if(iCascType == 1) Printf("Xi- current index = %n ", iCurrentLabelStack);
2062                         ncascperev++;
2063                 }// end if current particle = Xi-
2064
2065      }// This is the end of the loop on primaries
2066
2067      if(iCascType == 1) {
2068        fHistnXiMinusPerEv->Fill(ncascperev);
2069        fHistnXiMinusPerEvTot->Fill(ncascperevtot);
2070 //       Printf("N xi-tot = %n N xi-acc = %n\n", ncascperevtot, ncascperev);
2071      }
2072      if(iCascType == 2) {
2073        fHistnXiPlusPerEv->Fill(ncascperev);
2074        fHistnXiPlusPerEvTot->Fill(ncascperevtot);
2075      }
2076      if(iCascType == 3) {
2077        fHistnOmegaMinusPerEv->Fill(ncascperev);
2078        fHistnOmegaMinusPerEvTot->Fill(ncascperevtot);
2079      }
2080      if(iCascType == 4) {
2081        fHistnOmegaPlusPerEv->Fill(ncascperev);
2082        fHistnOmegaPlusPerEvTot->Fill(ncascperevtot);
2083      }
2084
2085
2086
2087
2088 // - Re-initialisation of the local TH1F pointers
2089 lHistEtaGenCasc         = 0x0;
2090 l3dHistGenPtVsGenYGenvsCent = 0x0;
2091 l3dHistGenPtVsGenYGenvsNtracks = 0x0;
2092
2093 lHistThetaGenCasc       = 0x0;
2094 l2dHistGenPtVsGenYFdbl  = 0x0;
2095 lHistThetaLambda        = 0x0;
2096 lHistThetaBach          = 0x0;
2097 lHistThetaBarDghter     = 0x0;
2098 lHistThetaMesDghter     = 0x0;
2099 lHistPtBach             = 0x0;
2100 lHistPtBarDghter        = 0x0;
2101 lHistPtMesDghter        = 0x0;
2102
2103 } // end of loop over the different types of cascades (Xi-+, Omega-+)
2104
2105
2106
2107 //__________________________________________________________________________
2108 // Part 3 - Loop over the reconstructed candidates
2109
2110 Int_t nAssoXiMinus = 0;
2111 Int_t nAssoXiPlus = 0;
2112 Int_t nAssoOmegaMinus = 0;
2113 Int_t nAssoOmegaPlus = 0;
2114
2115 for (Int_t iXi = 0; iXi < ncascades; iXi++) {// This is the begining of the Cascade loop
2116
2117         AliESDcascade *xiESD = lESDevent->GetCascade(iXi);
2118         if (!xiESD) continue;
2119
2120         // - Step II.1 : Connection to daughter tracks of the current cascade
2121         //-------------
2122
2123                 UInt_t lIdxPosXi        = (UInt_t) TMath::Abs( xiESD->GetPindex() );
2124                 UInt_t lIdxNegXi        = (UInt_t) TMath::Abs( xiESD->GetNindex() );
2125                 UInt_t lBachIdx         = (UInt_t) TMath::Abs( xiESD->GetBindex() );
2126                 // abs value not needed ; the index should always be positive (!= label ...)
2127
2128
2129         // FIXME : rejection of a double use of a daughter track (nothing but just a crosscheck of what is done in the cascade vertexer)
2130         if(lBachIdx == lIdxNegXi) {
2131                 AliWarning("Pb / Idx(Bach. track) = Idx(Neg. track) ... continue!"); continue;
2132         }
2133         if(lBachIdx == lIdxPosXi) {
2134                 AliWarning("Pb / Idx(Bach. track) = Idx(Pos. track) ... continue!"); continue;
2135         }
2136
2137         AliESDtrack *pTrackXi           = lESDevent->GetTrack( lIdxPosXi );
2138         AliESDtrack *nTrackXi           = lESDevent->GetTrack( lIdxNegXi );
2139         AliESDtrack *bachTrackXi        = lESDevent->GetTrack( lBachIdx  );
2140         if (!pTrackXi || !nTrackXi || !bachTrackXi ) {
2141                 Printf("ERROR: Could not retrieve one of the 3 daughter tracks of the cascade ...");
2142                 continue;
2143         }
2144
2145         Int_t lPosTPCClusters   = pTrackXi->GetTPCNcls();
2146         Int_t lNegTPCClusters   = nTrackXi->GetTPCNcls();
2147         Int_t lBachTPCClusters  = bachTrackXi->GetTPCNcls();
2148
2149                 // FIXME : rejection of a poor quality tracks
2150         if(fkQualityCutTPCrefit){
2151                 // 1 - Poor quality related to TPCrefit
2152                 ULong_t pStatus    = pTrackXi->GetStatus();
2153                 ULong_t nStatus    = nTrackXi->GetStatus();
2154                 ULong_t bachStatus = bachTrackXi->GetStatus();
2155                 if ((pStatus&AliESDtrack::kTPCrefit)    == 0) { AliWarning("Pb / V0 Pos. track has no TPCrefit ... continue!"); continue; }
2156                 if ((nStatus&AliESDtrack::kTPCrefit)    == 0) { AliWarning("Pb / V0 Neg. track has no TPCrefit ... continue!"); continue; }
2157                 if ((bachStatus&AliESDtrack::kTPCrefit) == 0) { AliWarning("Pb / Bach.   track has no TPCrefit ... continue!"); continue; }
2158         }
2159         if(fkQualityCut80TPCcls){
2160                 // 2 - Poor quality related to TPC clusters
2161                 if(lPosTPCClusters  < 80) { AliWarning("Pb / V0 Pos. track has less than 80 TPC clusters ... continue!"); continue; }
2162                 if(lNegTPCClusters  < 80) { AliWarning("Pb / V0 Neg. track has less than 80 TPC clusters ... continue!"); continue; }
2163                 if(lBachTPCClusters < 80) { AliWarning("Pb / Bach.   track has less than 80 TPC clusters ... continue!"); continue; }
2164         }
2165
2166         // - Step II.2 : Info over reconstructed cascades
2167         //-------------
2168
2169         Double_t lInvMassXiMinus    = 0.;
2170         Double_t lInvMassXiPlus     = 0.;
2171         Double_t lInvMassOmegaMinus = 0.;
2172         Double_t lInvMassOmegaPlus  = 0.;
2173
2174         Double_t lV0quality = 0.;
2175
2176         if( bachTrackXi->Charge() < 0 ) {
2177                 lV0quality = 0.;
2178                 xiESD->ChangeMassHypothesis(lV0quality , 3312);
2179                         // Calculate the effective mass of the Xi- candidate.
2180                         // pdg code 3312 = Xi-
2181                 lInvMassXiMinus = xiESD->GetEffMassXi();
2182
2183                 lV0quality = 0.;
2184                 xiESD->ChangeMassHypothesis(lV0quality , 3334);
2185                         // Calculate the effective mass of the Xi- candidate.
2186                         // pdg code 3334 = Omega-
2187                 lInvMassOmegaMinus = xiESD->GetEffMassXi();
2188
2189                 lV0quality = 0.;
2190                 xiESD->ChangeMassHypothesis(lV0quality , 3312);         // Back to default hyp.
2191
2192         }
2193
2194         if( bachTrackXi->Charge() >  0 ){
2195                 lV0quality = 0.;
2196                 xiESD->ChangeMassHypothesis(lV0quality , -3312);
2197                         // Calculate the effective mass of the Xi+ candidate.
2198                         // pdg code -3312 = Xi+
2199                 lInvMassXiPlus = xiESD->GetEffMassXi();
2200
2201                 lV0quality = 0.;
2202                 xiESD->ChangeMassHypothesis(lV0quality , -3334);
2203                         // Calculate the effective mass of the Xi+ candidate.
2204                         // pdg code -3334  = Omega+
2205                 lInvMassOmegaPlus = xiESD->GetEffMassXi();
2206
2207                 lV0quality = 0.;
2208                 xiESD->ChangeMassHypothesis(lV0quality , -3312);        // Back to "default" hyp.
2209         }
2210
2211         // Double_t lChi2Xi         = -1. ;
2212         Double_t lDcaXiDaughters            = -1. ;
2213         Double_t lDcaBachToPrimVertexXi     = -1. ;
2214         Double_t lXiCosineOfPointingAngle   = -1. ;
2215         Double_t lPosXi[3]                  = { -1000.0, -1000.0, -1000.0 };
2216         Double_t lXiRadius                  = -1000. ;
2217
2218         Double_t lInvMassLambdaAsCascDghter = 0.;
2219         Double_t lDcaV0DaughtersXi          = -1.;
2220         // Double_t lV0Chi2Xi               = -1. ;
2221         Double_t lV0CosineOfPointingAngleXi = -1.;
2222         Double_t lPosV0Xi[3]                = { -1000. , -1000., -1000. }; // Position of VO coming from cascade
2223         Double_t lV0RadiusXi                = -1000.;
2224         Double_t lDcaV0ToPrimVertexXi       = -1.;
2225         Double_t lDcaPosToPrimVertexXi      = -1.;
2226         Double_t lDcaNegToPrimVertexXi      = -1.;
2227
2228         Int_t    nTrackWithTPCrefitMultiplicity  =  0;
2229         Int_t    lSPDTrackletsMultiplicity       = -1;
2230
2231         //lChi2Xi                       = xiESD->GetChi2Xi();
2232         lDcaXiDaughters                 = xiESD->GetDcaXiDaughters();
2233         lDcaBachToPrimVertexXi          = TMath::Abs( bachTrackXi->GetD( lBestPrimaryVtxPos[0],
2234                                                                          lBestPrimaryVtxPos[1],
2235                                                                          lMagneticField  ) );
2236                                         // NOTE : AliExternalTrackParam::GetD returns an algebraic value
2237         lXiCosineOfPointingAngle        = xiESD->GetCascadeCosineOfPointingAngle( lBestPrimaryVtxPos[0],
2238                                                                                   lBestPrimaryVtxPos[1],
2239                                                                                   lBestPrimaryVtxPos[2] );
2240                                         // Take care : the best available vertex should be used (like in AliCascadeVertexer)
2241         xiESD->GetXYZcascade( lPosXi[0],  lPosXi[1], lPosXi[2] );
2242         lXiRadius                       = TMath::Sqrt( lPosXi[0]*lPosXi[0]  +  lPosXi[1]*lPosXi[1] );
2243         lInvMassLambdaAsCascDghter      = xiESD->GetEffMass();
2244                                         // This value shouldn't change, whatever the working hyp. is : Xi-, Xi+, Omega-, Omega+
2245         lDcaV0DaughtersXi               = xiESD->GetDcaV0Daughters();
2246         // lV0Chi2Xi                    = xiESD->GetChi2V0();
2247         lV0CosineOfPointingAngleXi      = xiESD->GetV0CosineOfPointingAngle( lBestPrimaryVtxPos[0],
2248                                                                              lBestPrimaryVtxPos[1],
2249                                                                              lBestPrimaryVtxPos[2] );
2250         xiESD->GetXYZ( lPosV0Xi[0],  lPosV0Xi[1], lPosV0Xi[2] );
2251         lV0RadiusXi                     = TMath::Sqrt( lPosV0Xi[0]*lPosV0Xi[0]  +  lPosV0Xi[1]*lPosV0Xi[1] );
2252
2253         lDcaV0ToPrimVertexXi            = xiESD->GetD( lBestPrimaryVtxPos[0],
2254                                                 lBestPrimaryVtxPos[1],
2255                                                 lBestPrimaryVtxPos[2] );
2256
2257         lDcaPosToPrimVertexXi           = TMath::Abs( pTrackXi   ->GetD( lBestPrimaryVtxPos[0],
2258                                                                          lBestPrimaryVtxPos[1],
2259                                                                          lMagneticField  )     );
2260
2261         lDcaNegToPrimVertexXi           = TMath::Abs( nTrackXi   ->GetD( lBestPrimaryVtxPos[0],
2262                                                                          lBestPrimaryVtxPos[1],
2263                                                                          lMagneticField  )     );
2264
2265                                // - II.Step 3' : extra-selection for cascade candidates
2266                 // Towards optimisation of AA selection
2267       if (fkExtraSelections) {
2268
2269         // if(lChi2Xi > 2000) continue; // in AliCascadeVertexer
2270         // if(lV0Chi2Xi > 2000) continue; // in AliV0vertexer
2271
2272         if (lDcaXiDaughters > 0.3) continue; // in AliCascadeVertexer
2273         if (lXiCosineOfPointingAngle < 0.999 ) continue; // in AliCascadeVertexer
2274         if (lDcaV0ToPrimVertexXi < 0.05) continue; // in AliCascadeVertexer
2275         if (lDcaBachToPrimVertexXi < 0.03) continue; // in AliCascadeVertexer
2276 ////      if (TMath::Abs(lInvMassLambdaAsCascDghter-1.11568) > 0.006 ) continue;  // in AliCascadeVertexer
2277
2278         if (lDcaV0DaughtersXi > 1.) continue; // in AliV0vertexer
2279         if (lV0CosineOfPointingAngleXi < 0.998) continue; // in AliV0vertexer
2280         if (lDcaPosToPrimVertexXi < 0.1) continue; // in AliV0vertexer
2281         if (lDcaNegToPrimVertexXi < 0.1) continue; // in AliV0vertexer
2282
2283
2284           if(lXiRadius < .9) continue; // in AliCascadeVertexer
2285 //        if(lXiRadius > 100) continue; // in AliCascadeVertexer
2286           if(lV0RadiusXi < 0.9) continue; // in AliV0vertexer
2287 //        if(lV0RadiusXi > 100) continue; // in AliV0vertexer
2288
2289       }
2290
2291         Double_t lChargeXi = xiESD->Charge();
2292
2293         if( lChargeXi < 0 )     fHistMassXiMinus        ->Fill( lInvMassXiMinus );
2294         if( lChargeXi > 0 )     fHistMassXiPlus         ->Fill( lInvMassXiPlus );
2295         if( lChargeXi < 0 )     fHistMassOmegaMinus     ->Fill( lInvMassOmegaMinus );
2296         if( lChargeXi > 0 )     fHistMassOmegaPlus      ->Fill( lInvMassOmegaPlus );
2297
2298         // - Step II.3 : PID info
2299
2300         //-------------
2301
2302
2303         // 3.1 - PID Information
2304
2305         Bool_t   lIsPosInXiProton      = kFALSE;
2306         Bool_t   lIsPosInXiPion        = kFALSE;
2307         Bool_t   lIsPosInOmegaProton   = kFALSE;
2308         Bool_t   lIsPosInOmegaPion     = kFALSE;
2309
2310         Bool_t   lIsNegInXiProton      = kFALSE;
2311         Bool_t   lIsNegInXiPion        = kFALSE;
2312         Bool_t   lIsNegInOmegaProton   = kFALSE;
2313         Bool_t   lIsNegInOmegaPion     = kFALSE;
2314
2315         Bool_t   lIsBachelorKaon       = kFALSE;
2316         Bool_t   lIsBachelorPion       = kFALSE;
2317
2318         Bool_t   lIsBachelorKaonForTPC = kFALSE; // For ESD only ...//FIXME : wait for availability in AOD
2319         Bool_t   lIsBachelorPionForTPC = kFALSE; // For ESD only ...
2320         Bool_t   lIsNegPionForTPC      = kFALSE; // For ESD only ...
2321         Bool_t   lIsPosPionForTPC      = kFALSE; // For ESD only ...
2322         Bool_t   lIsNegProtonForTPC    = kFALSE; // For ESD only ...
2323         Bool_t   lIsPosProtonForTPC    = kFALSE; // For ESD only ...
2324
2325         // 3.1.A - Combined PID
2326         // Reasonable guess for the priors for the cascade track sample (e-, mu, pi, K, p)
2327         Double_t lPriorsGuessXi[5]    = {0, 0, 2, 0, 1};
2328         Double_t lPriorsGuessOmega[5] = {0, 0, 1, 1, 1};
2329
2330         // Combined VO-positive-daughter PID
2331         AliPID pPidXi;          pPidXi.SetPriors(    lPriorsGuessXi    );
2332         AliPID pPidOmega;       pPidOmega.SetPriors( lPriorsGuessOmega );
2333
2334         if( pTrackXi->IsOn(AliESDtrack::kESDpid) ){  // Combined PID exists
2335                 Double_t r[10] = {0.}; pTrackXi->GetESDpid(r);
2336                 pPidXi.SetProbabilities(r);
2337                 pPidOmega.SetProbabilities(r);
2338
2339                 // Check if the V0 positive track is a proton (case for Xi-)
2340                 Double_t pproton = pPidXi.GetProbability(AliPID::kProton);
2341                 if (pproton > pPidXi.GetProbability(AliPID::kElectron) &&
2342                     pproton > pPidXi.GetProbability(AliPID::kMuon)     &&
2343                     pproton > pPidXi.GetProbability(AliPID::kPion)     &&
2344                     pproton > pPidXi.GetProbability(AliPID::kKaon)     )     lIsPosInXiProton = kTRUE;
2345
2346                 // Check if the V0 positive track is a pi+ (case for Xi+)
2347                 Double_t ppion = pPidXi.GetProbability(AliPID::kPion);
2348                 if (ppion > pPidXi.GetProbability(AliPID::kElectron) &&
2349                     ppion > pPidXi.GetProbability(AliPID::kMuon)     &&
2350                     ppion > pPidXi.GetProbability(AliPID::kKaon)     &&
2351                     ppion > pPidXi.GetProbability(AliPID::kProton)   )     lIsPosInXiPion = kTRUE;
2352
2353
2354                 // Check if the V0 positive track is a proton (case for Omega-)
2355                 pproton = 0.;
2356                     pproton = pPidOmega.GetProbability(AliPID::kProton);
2357                 if (pproton > pPidOmega.GetProbability(AliPID::kElectron) &&
2358                     pproton > pPidOmega.GetProbability(AliPID::kMuon)     &&
2359                     pproton > pPidOmega.GetProbability(AliPID::kPion)     &&
2360                     pproton > pPidOmega.GetProbability(AliPID::kKaon)     )  lIsPosInOmegaProton = kTRUE;
2361
2362                 // Check if the V0 positive track is a pi+ (case for Omega+)
2363                 ppion = 0.;
2364                     ppion = pPidOmega.GetProbability(AliPID::kPion);
2365                 if (ppion > pPidOmega.GetProbability(AliPID::kElectron) &&
2366                     ppion > pPidOmega.GetProbability(AliPID::kMuon)     &&
2367                     ppion > pPidOmega.GetProbability(AliPID::kKaon)     &&
2368                     ppion > pPidOmega.GetProbability(AliPID::kProton)   )    lIsPosInOmegaPion = kTRUE;
2369
2370         }// end if V0 positive track with existing combined PID
2371
2372
2373         // Combined VO-negative-daughter PID
2374         AliPID nPidXi;          nPidXi.SetPriors(    lPriorsGuessXi    );
2375         AliPID nPidOmega;       nPidOmega.SetPriors( lPriorsGuessOmega );
2376
2377         if( nTrackXi->IsOn(AliESDtrack::kESDpid) ){  // Combined PID exists
2378                 Double_t r[10] = {0.}; nTrackXi->GetESDpid(r);
2379                 nPidXi.SetProbabilities(r);
2380                 nPidOmega.SetProbabilities(r);
2381
2382                 // Check if the V0 negative track is a pi- (case for Xi-)
2383                 Double_t ppion = nPidXi.GetProbability(AliPID::kPion);
2384                 if (ppion > nPidXi.GetProbability(AliPID::kElectron) &&
2385                     ppion > nPidXi.GetProbability(AliPID::kMuon)     &&
2386                     ppion > nPidXi.GetProbability(AliPID::kKaon)     &&
2387                     ppion > nPidXi.GetProbability(AliPID::kProton)   )     lIsNegInXiPion = kTRUE;
2388
2389                 // Check if the V0 negative track is an anti-proton (case for Xi+)
2390                 Double_t pproton = nPidXi.GetProbability(AliPID::kProton);
2391                 if (pproton > nPidXi.GetProbability(AliPID::kElectron) &&
2392                     pproton > nPidXi.GetProbability(AliPID::kMuon)     &&
2393                     pproton > nPidXi.GetProbability(AliPID::kPion)     &&
2394                     pproton > nPidXi.GetProbability(AliPID::kKaon)     )     lIsNegInXiProton = kTRUE;
2395
2396                 // Check if the V0 negative track is a pi- (case for Omega-)
2397                 ppion = 0.;
2398                     ppion = nPidOmega.GetProbability(AliPID::kPion);
2399                 if (ppion > nPidOmega.GetProbability(AliPID::kElectron) &&
2400                     ppion > nPidOmega.GetProbability(AliPID::kMuon)     &&
2401                     ppion > nPidOmega.GetProbability(AliPID::kKaon)     &&
2402                     ppion > nPidOmega.GetProbability(AliPID::kProton)   )    lIsNegInOmegaPion = kTRUE;
2403
2404                 // Check if the V0 negative track is an anti-proton (case for Omega+)
2405                 pproton = 0.;
2406                     pproton = nPidOmega.GetProbability(AliPID::kProton);
2407                 if (pproton > nPidOmega.GetProbability(AliPID::kElectron) &&
2408                     pproton > nPidOmega.GetProbability(AliPID::kMuon)     &&
2409                     pproton > nPidOmega.GetProbability(AliPID::kPion)     &&
2410                     pproton > nPidOmega.GetProbability(AliPID::kKaon)     )  lIsNegInOmegaProton = kTRUE;
2411
2412         }// end if V0 negative track with existing combined PID
2413
2414
2415         // Combined bachelor PID
2416         AliPID bachPidXi;       bachPidXi.SetPriors(    lPriorsGuessXi    );
2417         AliPID bachPidOmega;    bachPidOmega.SetPriors( lPriorsGuessOmega );
2418
2419         Double_t ppionBach = 0.0, pkaonBach = 0.0;
2420
2421         if( bachTrackXi->IsOn(AliESDtrack::kESDpid) ){  // Combined PID exists
2422                 Double_t r[10] = {0.}; bachTrackXi->GetESDpid(r);
2423                 bachPidXi.SetProbabilities(r);
2424                 bachPidOmega.SetProbabilities(r);
2425                 // Check if the bachelor track is a pion
2426                     ppionBach = bachPidXi.GetProbability(AliPID::kPion);
2427                 if (ppionBach > bachPidXi.GetProbability(AliPID::kElectron) &&
2428                     ppionBach > bachPidXi.GetProbability(AliPID::kMuon)     &&
2429                     ppionBach > bachPidXi.GetProbability(AliPID::kKaon)     &&
2430                     ppionBach > bachPidXi.GetProbability(AliPID::kProton)   )     lIsBachelorPion = kTRUE;
2431                 // Check if the bachelor track is a kaon
2432                     pkaonBach = bachPidOmega.GetProbability(AliPID::kKaon);
2433                 if (pkaonBach > bachPidOmega.GetProbability(AliPID::kElectron) &&
2434                     pkaonBach > bachPidOmega.GetProbability(AliPID::kMuon)     &&
2435                     pkaonBach > bachPidOmega.GetProbability(AliPID::kPion)     &&
2436                     pkaonBach > bachPidOmega.GetProbability(AliPID::kProton)   )  lIsBachelorKaon = kTRUE;
2437         }// end if bachelor track with existing combined PID
2438
2439
2440         // 3.1.B - TPC PID : 4-sigma bands on Bethe-Bloch curve
2441
2442         // Here for the new PID object I put fPIDResponse instead of fESDpid
2443         // Bachelor
2444         if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( bachTrackXi,AliPID::kKaon)) < 4) lIsBachelorKaonForTPC = kTRUE;
2445         if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( bachTrackXi,AliPID::kPion)) < 4) lIsBachelorPionForTPC = kTRUE;
2446
2447         // Negative V0 daughter
2448         if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( nTrackXi,AliPID::kPion   )) < 4) lIsNegPionForTPC   = kTRUE;
2449         if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( nTrackXi,AliPID::kProton )) < 4) lIsNegProtonForTPC = kTRUE;
2450
2451         // Positive V0 daughter
2452         if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( pTrackXi,AliPID::kPion   )) < 4) lIsPosPionForTPC   = kTRUE;
2453         if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( pTrackXi,AliPID::kProton )) < 4) lIsPosProtonForTPC = kTRUE;
2454
2455 /*
2456         const AliExternalTrackParam *pInnerWallTrackXi    = pTrackXi    ->GetInnerParam(); // Do not use GetTPCInnerWall
2457         const AliExternalTrackParam *nInnerWallTrackXi    = nTrackXi    ->GetInnerParam();
2458         const AliExternalTrackParam *bachInnerWallTrackXi = bachTrackXi ->GetInnerParam();
2459         if(pInnerWallTrackXi && nInnerWallTrackXi && bachInnerWallTrackXi ){
2460
2461                 Double_t pMomInnerWall    = pInnerWallTrackXi   ->GetP();
2462                 Double_t nMomInnerWall    = nInnerWallTrackXi   ->GetP();
2463                 Double_t bachMomInnerWall = bachInnerWallTrackXi->GetP();
2464
2465                 // Bachelor
2466                 if (TMath::Abs(fESDpid->NumberOfSigmasTPC( bachTrackXi,AliPID::kPion)) < 3)                              lIsBachelorPionForTPC = kTRUE;
2467                 if (bachMomInnerWall < 0.350  && TMath::Abs(fESDpid->NumberOfSigmasTPC( bachTrackXi,AliPID::kKaon)) < 5) lIsBachelorKaonForTPC = kTRUE;
2468                 if (bachMomInnerWall > 0.350  && TMath::Abs(fESDpid->NumberOfSigmasTPC( bachTrackXi,AliPID::kKaon)) < 3) lIsBachelorKaonForTPC = kTRUE;
2469
2470                 // Negative V0 daughter
2471                 if (TMath::Abs(fESDpid->NumberOfSigmasTPC( nTrackXi,AliPID::kPion   )) < 3  )                           lIsNegPionForTPC   = kTRUE;
2472                 if (nMomInnerWall < 0.6  && TMath::Abs(fESDpid->NumberOfSigmasTPC( nTrackXi,AliPID::kProton ) ) < 5 )   lIsNegProtonForTPC = kTRUE;
2473                 if (nMomInnerWall > 0.6  && TMath::Abs(fESDpid->NumberOfSigmasTPC( nTrackXi,AliPID::kProton ) ) < 3 )   lIsNegProtonForTPC = kTRUE;
2474
2475                 // Positive V0 daughter
2476                 if (TMath::Abs(fESDpid->NumberOfSigmasTPC( pTrackXi,AliPID::kPion   )) < 3 )                            lIsPosPionForTPC   = kTRUE;
2477                 if (pMomInnerWall < 0.6  && TMath::Abs(fESDpid->NumberOfSigmasTPC( pTrackXi,AliPID::kProton )) < 5)     lIsPosProtonForTPC = kTRUE;
2478                 if (pMomInnerWall > 0.6  && TMath::Abs(fESDpid->NumberOfSigmasTPC( pTrackXi,AliPID::kProton )) < 3)     lIsPosProtonForTPC = kTRUE;
2479
2480         }
2481 */
2482
2483         // Combined PID TH1s
2484         if( lChargeXi < 0 && lIsBachelorPion )    fHistMassWithCombPIDXiMinus     ->Fill( lInvMassXiMinus    );
2485         if( lChargeXi > 0 && lIsBachelorPion )    fHistMassWithCombPIDXiPlus      ->Fill( lInvMassXiPlus     );
2486         if( lChargeXi < 0 && lIsBachelorKaon )    fHistMassWithCombPIDOmegaMinus  ->Fill( lInvMassOmegaMinus );
2487         if( lChargeXi > 0 && lIsBachelorKaon )    fHistMassWithCombPIDOmegaPlus   ->Fill( lInvMassOmegaPlus  );
2488
2489
2490         // 3.2 - PID proba Vs Pt(Bach)
2491         Int_t      lblBachForPID  = (Int_t) TMath::Abs( bachTrackXi->GetLabel() );
2492         TParticle* mcBachForPID   = lMCstack->Particle( lblBachForPID );
2493         Double_t   lmcPtBach      = mcBachForPID->Pt();
2494
2495         if(lIsBachelorPion)   f2dHistPIDprobaPionVsMCPtBach->Fill( lmcPtBach, ppionBach );
2496         if(lIsBachelorKaon)   f2dHistPIDprobaKaonVsMCPtBach->Fill( lmcPtBach, pkaonBach );
2497
2498
2499         // 3.3 - MC perfect PID
2500         Bool_t   lIsBachelorMCPiMinus  = kFALSE;
2501         Bool_t   lIsBachelorMCPiPlus   = kFALSE;
2502         Bool_t   lIsBachelorMCKMinus   = kFALSE;
2503         Bool_t   lIsBachelorMCKPlus    = kFALSE;
2504
2505         if( mcBachForPID->GetPdgCode() == -211) lIsBachelorMCPiMinus = kTRUE;
2506         if( mcBachForPID->GetPdgCode() ==  211) lIsBachelorMCPiPlus  = kTRUE;
2507         if( mcBachForPID->GetPdgCode() == -321) lIsBachelorMCKMinus  = kTRUE;
2508         if( mcBachForPID->GetPdgCode() ==  321) lIsBachelorMCKPlus   = kTRUE;
2509
2510         if( lChargeXi < 0 && lIsBachelorMCPiMinus )    fHistMassWithMcPIDXiMinus     ->Fill( lInvMassXiMinus );
2511         if( lChargeXi > 0 && lIsBachelorMCPiPlus  )    fHistMassWithMcPIDXiPlus      ->Fill( lInvMassXiPlus );
2512         if( lChargeXi < 0 && lIsBachelorMCKMinus  )    fHistMassWithMcPIDOmegaMinus  ->Fill( lInvMassOmegaMinus );
2513         if( lChargeXi > 0 && lIsBachelorMCKPlus   )    fHistMassWithMcPIDOmegaPlus   ->Fill( lInvMassOmegaPlus );
2514
2515
2516         // - Step II.4 : MC association (care : lots of "continue;" below this line)
2517         //-------------
2518
2519         Bool_t lAssoXiMinus    = kFALSE;
2520         Bool_t lAssoXiPlus     = kFALSE;
2521         Bool_t lAssoOmegaMinus = kFALSE;
2522         Bool_t lAssoOmegaPlus  = kFALSE;
2523
2524
2525         if(fDebug > 5)
2526                 cout    << "MC EventNumber : " << lMCevent->Header()->GetEvent()
2527                         << " / MC event Number in Run : " << lMCevent->Header()->GetEventNrInRun() << endl;
2528
2529         // - Step 4.1 : level of the V0 daughters
2530
2531         Int_t lblPosV0Dghter = (Int_t) TMath::Abs( pTrackXi->GetLabel() );
2532                 // Abs value = needed ! question of quality track association ... (negative label when at least one cluster in the track is from a different particle)
2533         Int_t lblNegV0Dghter = (Int_t) TMath::Abs( nTrackXi->GetLabel() );
2534         TParticle* mcPosV0Dghter = lMCstack->Particle( lblPosV0Dghter );
2535         TParticle* mcNegV0Dghter = lMCstack->Particle( lblNegV0Dghter );
2536
2537
2538         // - Step 4.2 : level of the Xi daughters
2539
2540         Int_t lblMotherPosV0Dghter = mcPosV0Dghter->GetFirstMother() ;
2541         Int_t lblMotherNegV0Dghter = mcNegV0Dghter->GetFirstMother();
2542
2543                 if( lblMotherPosV0Dghter != lblMotherNegV0Dghter) continue; // same mother
2544                 if( lblMotherPosV0Dghter < 0 ) continue; // this particle is primary, no mother
2545                 if( lblMotherNegV0Dghter < 0 ) continue;
2546
2547
2548                 // mothers = Lambda candidate ... a priori
2549
2550         TParticle* mcMotherPosV0Dghter = lMCstack->Particle( lblMotherPosV0Dghter );
2551         TParticle* mcMotherNegV0Dghter = lMCstack->Particle( lblMotherNegV0Dghter );  // MN: redundant?? already checked that labels are the same...-->same part from stack
2552
2553         Int_t      lblBach  = (Int_t) TMath::Abs( bachTrackXi->GetLabel() );
2554         TParticle* mcBach   = lMCstack->Particle( lblBach );
2555
2556
2557         // - Step 4.3 : level of Xi candidate
2558
2559         Int_t lblGdMotherPosV0Dghter =   mcMotherPosV0Dghter->GetFirstMother() ;
2560         Int_t lblGdMotherNegV0Dghter =   mcMotherNegV0Dghter->GetFirstMother() ;
2561
2562                 if( lblGdMotherPosV0Dghter != lblGdMotherNegV0Dghter ) continue;
2563                 if( lblGdMotherPosV0Dghter < 0 ) continue; // primary lambda ...
2564                 if( lblGdMotherNegV0Dghter < 0 ) continue; // primary lambda ...
2565
2566                 // Gd mothers = Xi candidate ... a priori
2567
2568         TParticle* mcGdMotherPosV0Dghter = lMCstack->Particle( lblGdMotherPosV0Dghter );
2569         TParticle* mcGdMotherNegV0Dghter = lMCstack->Particle( lblGdMotherNegV0Dghter );
2570
2571         Int_t lblMotherBach = (Int_t) TMath::Abs( mcBach->GetFirstMother()  );
2572
2573         if( lblMotherBach != lblGdMotherPosV0Dghter ) continue; //same mother for bach and V0 daughters
2574
2575         TParticle* mcMotherBach = lMCstack->Particle( lblMotherBach );
2576
2577         // Check if cascade is primary
2578
2579         if (!(lMCstack->IsPhysicalPrimary(lblMotherBach))) continue;
2580
2581         // - Step 4.4 : Manage boolean for association
2582
2583         if( mcMotherBach                ->GetPdgCode() ==   3312 &&
2584             mcGdMotherPosV0Dghter       ->GetPdgCode() ==   3312 &&
2585             mcGdMotherNegV0Dghter       ->GetPdgCode() ==   3312)    {  lAssoXiMinus = kTRUE;
2586                                                                       //  Printf("Xi- asso current index = %n ", lblGdMotherPosV0Dghter);
2587                                                                         nAssoXiMinus++; }
2588
2589         else if( mcMotherBach           ->GetPdgCode() ==  -3312 &&
2590             mcGdMotherPosV0Dghter       ->GetPdgCode() ==  -3312 &&
2591             mcGdMotherNegV0Dghter       ->GetPdgCode() ==  -3312)    {  lAssoXiPlus = kTRUE;
2592                                                                         nAssoXiPlus++; }
2593
2594         else if( mcMotherBach           ->GetPdgCode() ==   3334 &&
2595             mcGdMotherPosV0Dghter       ->GetPdgCode() ==   3334 &&
2596             mcGdMotherNegV0Dghter       ->GetPdgCode() ==   3334)    {  lAssoOmegaMinus = kTRUE;
2597                                                                         nAssoOmegaMinus++; }
2598
2599         else if( mcMotherBach           ->GetPdgCode() ==  -3334 &&
2600             mcGdMotherPosV0Dghter       ->GetPdgCode() ==  -3334 &&
2601             mcGdMotherNegV0Dghter       ->GetPdgCode() ==  -3334)    {  lAssoOmegaPlus = kTRUE;
2602                                                                         nAssoOmegaPlus++; }
2603
2604
2605
2606         if(!lAssoXiMinus && !lAssoXiPlus && !lAssoOmegaMinus && !lAssoOmegaPlus) continue; // no association, skip the rest of the code
2607         // If a proper association  exists ...
2608
2609         if(fDebug > 4){
2610                 cout << "XiMinus    = " << lAssoXiMinus    << endl;
2611                 cout << "XiPlus     = " << lAssoXiPlus     << endl;
2612                 cout << "OmegaMinus = " << lAssoOmegaMinus << endl;
2613                 cout << "OmegaPlus  = " << lAssoOmegaPlus  << endl
2614                      << "----"          << endl;
2615         }
2616
2617
2618         if(fDebug > 5){
2619                 cout << endl;
2620                 cout << "- V0 daughters - " << endl;
2621                 cout << "     + V0 Pos. / Label : " << lblPosV0Dghter
2622                 << " - Pdg Code : " << mcPosV0Dghter->GetTitle() << endl;
2623                 cout << "     - V0 Neg. / Label : " << lblNegV0Dghter
2624                 << " - Pdg Code : " << mcNegV0Dghter->GetTitle() << endl;
2625
2626                 cout << "- Xi daughters - " << endl;
2627                 cout << "     + V0 Pos. mother / Label : " << lblMotherPosV0Dghter
2628                 << " - Pdg Code : " << mcMotherPosV0Dghter->GetTitle() << endl;
2629                 cout << "     - V0 Neg. mother / Label : " << lblMotherNegV0Dghter
2630                 << " - Pdg Code : " << mcMotherNegV0Dghter->GetTitle() << endl;
2631
2632                 cout << "     --  Bach. / Label :" << lblBach
2633                 << " -  Pdg Code : " << mcBach->GetTitle() << endl;
2634
2635                 cout << "- Xi candidate -" << endl;
2636                 cout << "    +  V0 Pos. Gd Mother / Label : " << lblGdMotherPosV0Dghter
2637                 << " - Pdg Code : " << mcGdMotherPosV0Dghter->GetTitle() << endl;
2638                 cout << "    -  V0 Neg. Gd Mother / Label : "  << lblGdMotherNegV0Dghter
2639                 << " - Pdg Code : "<< mcGdMotherNegV0Dghter->GetTitle() << endl;
2640
2641                 cout << "    --  Mother Bach. / Label : " << lblMotherBach
2642                 << " - Pdg Code    : " << mcMotherBach->GetTitle() << endl;
2643                 cout << endl;
2644         }
2645
2646
2647         // - Step 5 : Plots around the cascade candidates associated with MC
2648         //-------------
2649
2650         Double_t lmcPt             = mcMotherBach->Pt();
2651         Double_t lmcRapCasc        = 0.5*TMath::Log( (mcMotherBach->Energy() + mcMotherBach->Pz()) /
2652                                                      (mcMotherBach->Energy() - mcMotherBach->Pz() +1.e-13) );
2653         Double_t lmcEta            = mcMotherBach->Eta();
2654         Double_t lmcTransvRadius   = mcBach->R(); // to get the decay point of Xi, = the production vertex of Bachelor ...
2655
2656         TVector3 lmcTVect3Mom( mcMotherBach->Px(), mcMotherBach->Py(), mcMotherBach->Pz() );
2657
2658         Double_t lrecoPt           = xiESD->Pt();
2659         Double_t lrecoTransvRadius = TMath::Sqrt( xiESD->Xv() * xiESD->Xv() + xiESD->Yv() * xiESD->Yv() );
2660
2661         TVector3 lrecoTVect3Mom( xiESD->Px(), xiESD->Py(), xiESD->Pz()  );
2662         Double_t lDeltaPhiMcReco   = lmcTVect3Mom.DeltaPhi( lrecoTVect3Mom ) * 180.0/TMath::Pi();
2663
2664
2665         // - Histos for the cascade candidates associated with MC
2666
2667         if( lChargeXi < 0 && lAssoXiMinus){
2668                 fHistAsMCMassXiMinus          ->Fill( lInvMassXiMinus  );
2669                 if(lIsBachelorPion)     f2dHistAsMCandCombPIDGenPtVsGenYXiMinus->Fill( lmcPt, lmcRapCasc );
2670                 f2dHistAsMCGenPtVsGenYXiMinus ->Fill( lmcPt, lmcRapCasc);
2671                 fHistAsMCGenEtaXiMinus        ->Fill( lmcEta           );
2672                 f2dHistAsMCResPtXiMinus       ->Fill( lmcPt,           (lrecoPt - lmcPt)/ lmcPt );
2673                 f2dHistAsMCResRXiMinus        ->Fill( lmcTransvRadius, (lrecoTransvRadius - lmcTransvRadius)/ lmcTransvRadius    );
2674                 f2dHistAsMCResPhiXiMinus      ->Fill( lmcPt, lDeltaPhiMcReco );
2675         }
2676
2677         else if( lChargeXi > 0 && lAssoXiPlus){
2678                 fHistAsMCMassXiPlus           ->Fill( lInvMassXiPlus   );
2679                 if(lIsBachelorPion)     f2dHistAsMCandCombPIDGenPtVsGenYXiPlus->Fill( lmcPt, lmcRapCasc );
2680                 f2dHistAsMCGenPtVsGenYXiPlus  ->Fill( lmcPt, lmcRapCasc);
2681                 fHistAsMCGenEtaXiPlus         ->Fill( lmcEta           );
2682                 f2dHistAsMCResPtXiPlus        ->Fill( lmcPt,           (lrecoPt - lmcPt)/ lmcPt );
2683                 f2dHistAsMCResRXiPlus         ->Fill( lmcTransvRadius, (lrecoTransvRadius - lmcTransvRadius)/ lmcTransvRadius    );
2684                 f2dHistAsMCResPhiXiPlus       ->Fill( lmcPt, lDeltaPhiMcReco );
2685         }
2686
2687         else if( lChargeXi < 0 && lAssoOmegaMinus){
2688                 fHistAsMCMassOmegaMinus          ->Fill( lInvMassOmegaMinus );
2689                 if(lIsBachelorKaon)     f2dHistAsMCandCombPIDGenPtVsGenYOmegaMinus->Fill( lmcPt, lmcRapCasc );
2690                 f2dHistAsMCGenPtVsGenYOmegaMinus ->Fill( lmcPt, lmcRapCasc  );
2691                 fHistAsMCGenEtaOmegaMinus        ->Fill( lmcEta             );
2692                 f2dHistAsMCResPtOmegaMinus       ->Fill( lmcPt,           (lrecoPt - lmcPt)/ lmcPt );
2693                 f2dHistAsMCResROmegaMinus        ->Fill( lmcTransvRadius, (lrecoTransvRadius - lmcTransvRadius)/ lmcTransvRadius    );
2694                 f2dHistAsMCResPhiOmegaMinus      ->Fill( lmcPt, lDeltaPhiMcReco );
2695         }
2696
2697         else if( lChargeXi > 0 && lAssoOmegaPlus){
2698                 fHistAsMCMassOmegaPlus           ->Fill( lInvMassOmegaPlus );
2699                 if(lIsBachelorKaon)     f2dHistAsMCandCombPIDGenPtVsGenYOmegaPlus->Fill( lmcPt, lmcRapCasc );
2700                 f2dHistAsMCGenPtVsGenYOmegaPlus  ->Fill( lmcPt, lmcRapCasc   );
2701                 fHistAsMCGenEtaOmegaPlus         ->Fill( lmcEta            );
2702                 f2dHistAsMCResPtOmegaPlus        ->Fill( lmcPt,           (lrecoPt - lmcPt)/ lmcPt );
2703                 f2dHistAsMCResROmegaPlus         ->Fill( lmcTransvRadius, (lrecoTransvRadius - lmcTransvRadius)/ lmcTransvRadius    );
2704                 f2dHistAsMCResPhiOmegaPlus       ->Fill( lmcPt, lDeltaPhiMcReco );
2705         }
2706
2707
2708         // - Step 6 : Containers = Cascade cuts + PID
2709         //-------------
2710       if (fUseCFCont) {
2711
2712         nTrackWithTPCrefitMultiplicity  = DoESDTrackWithTPCrefitMultiplicity(lESDevent);  // FIXME : variable which is not used anymore at the moment ...
2713                                                                                           //    ->  keep it while the task is still under development.
2714
2715
2716         const AliMultiplicity *lAliMult = lESDevent->GetMultiplicity();
2717         if(fAnalysisType == "ESD") lSPDTrackletsMultiplicity       = lAliMult->GetNumberOfTracklets();
2718         else if(fAnalysisType == "AOD") lSPDTrackletsMultiplicity = lAODevent->GetTracklets()->GetNumberOfTracklets();
2719
2720
2721         // 6.3 - Filling the AliCFContainer (optimisation of topological selections + systematics)
2722         Double_t lContainerCutVars[19] = {0.0};
2723
2724         lContainerCutVars[0]  = lDcaXiDaughters;
2725         lContainerCutVars[1]  = lDcaBachToPrimVertexXi;
2726         lContainerCutVars[2]  = lXiCosineOfPointingAngle;
2727         lContainerCutVars[3]  = lXiRadius;
2728         lContainerCutVars[4]  = lInvMassLambdaAsCascDghter;
2729         lContainerCutVars[5]  = lDcaV0DaughtersXi;
2730         lContainerCutVars[6]  = lV0CosineOfPointingAngleXi;
2731         lContainerCutVars[7]  = lV0RadiusXi;
2732         lContainerCutVars[8]  = lDcaV0ToPrimVertexXi;
2733         lContainerCutVars[9]  = lDcaPosToPrimVertexXi;
2734         lContainerCutVars[10] = lDcaNegToPrimVertexXi;
2735
2736         lContainerCutVars[13] = lmcPt;
2737
2738         lContainerCutVars[16] = lBestPrimaryVtxPos[2];
2739         if (fCollidingSystems) lContainerCutVars[17] = lcentrality;
2740         else lContainerCutVars[17] = lSPDTrackletsMultiplicity;
2741         lContainerCutVars[18] = nTrackPrimaryMultiplicity;
2742
2743         // All cases should be covered below
2744         if( lChargeXi < 0 && lAssoXiMinus    ) {
2745                 lContainerCutVars[11] = lInvMassXiMinus;
2746                 lContainerCutVars[12] = lInvMassOmegaMinus;//1.63;
2747                 lContainerCutVars[14] = lmcRapCasc;
2748                 lContainerCutVars[15] = -1.;
2749                 if ( lIsBachelorPionForTPC   && lIsPosProtonForTPC    && lIsNegPionForTPC )
2750                   fCFContAsCascadeCuts->Fill(lContainerCutVars,0); // for Xi-
2751         }
2752         if( lChargeXi > 0 && lAssoXiPlus     ){
2753                 lContainerCutVars[11] = lInvMassXiPlus;
2754                 lContainerCutVars[12] = lInvMassOmegaPlus;//1.26;
2755                 lContainerCutVars[14] = lmcRapCasc;
2756                 lContainerCutVars[15] = -1.;
2757                 if ( lIsBachelorPionForTPC   && lIsNegProtonForTPC    && lIsPosPionForTPC )
2758                   fCFContAsCascadeCuts->Fill(lContainerCutVars,1); // for Xi+
2759         }
2760         if( lChargeXi < 0 && lAssoOmegaMinus )  {
2761                 lContainerCutVars[11] = lInvMassXiMinus;//1.63;
2762                 lContainerCutVars[12] = lInvMassOmegaMinus;
2763                 lContainerCutVars[14] = -1.;
2764                 lContainerCutVars[15] = lmcRapCasc;
2765                 if ( lIsBachelorKaonForTPC   && lIsPosProtonForTPC    && lIsNegPionForTPC )
2766                   fCFContAsCascadeCuts->Fill(lContainerCutVars,2); // for Omega-
2767         }
2768         if( lChargeXi > 0 && lAssoOmegaPlus  ){
2769                 lContainerCutVars[11] = lInvMassXiPlus;//1.26;
2770                 lContainerCutVars[12] = lInvMassOmegaPlus;
2771                 lContainerCutVars[14] = -1.;
2772                 lContainerCutVars[15] = lmcRapCasc;
2773                 if ( lIsBachelorKaonForTPC   && lIsNegProtonForTPC    && lIsPosPionForTPC )
2774                   fCFContAsCascadeCuts->Fill(lContainerCutVars,3); // for Omega+
2775         }
2776
2777
2778         // 6.4 - Filling the AliCFContainers related to PID
2779
2780         Double_t lContainerPIDVars[4] = {0.0};
2781
2782
2783         // Xi Minus
2784         if( lChargeXi < 0 && lAssoXiMinus ) {
2785                 lContainerPIDVars[0] = lmcPt              ;
2786                 lContainerPIDVars[1] = lInvMassXiMinus    ;
2787                 lContainerPIDVars[2] = lmcRapCasc         ;
2788                 if (fCollidingSystems) lContainerPIDVars[3] = lcentrality;
2789                 else lContainerPIDVars[3] = lSPDTrackletsMultiplicity ;
2790
2791                 // No PID
2792                         fCFContCascadePIDAsXiMinus->Fill(lContainerPIDVars, 0); // No PID
2793                 // TPC PID
2794                 if( lIsBachelorPionForTPC  )
2795                         fCFContCascadePIDAsXiMinus->Fill(lContainerPIDVars, 1); // TPC PID / 4-#sigma cut on Bachelor track
2796
2797                 if( lIsBachelorPionForTPC &&
2798                     lIsPosProtonForTPC     )
2799                         fCFContCascadePIDAsXiMinus->Fill(lContainerPIDVars, 2); // TPC PID / 4-#sigma cut on Bachelor+Baryon tracks
2800
2801                 if( lIsBachelorPionForTPC &&
2802                     lIsPosProtonForTPC    &&
2803                     lIsNegPionForTPC       )
2804                         fCFContCascadePIDAsXiMinus->Fill(lContainerPIDVars, 3); // TPC PID / 4-#sigma cut on Bachelor+Baryon+Meson tracks
2805
2806                 // Combined PID
2807                 if( lIsBachelorPion        )
2808                         fCFContCascadePIDAsXiMinus->Fill(lContainerPIDVars, 4); // Comb. PID / Bachelor
2809
2810                 if( lIsBachelorPion       &&
2811                     lIsPosInXiProton    )
2812                         fCFContCascadePIDAsXiMinus->Fill(lContainerPIDVars, 5); // Comb. PID / Bachelor+Baryon
2813
2814                 if(lIsBachelorPion     &&
2815                    lIsPosInXiProton &&
2816                    lIsNegInXiPion    )
2817                         fCFContCascadePIDAsXiMinus->Fill(lContainerPIDVars, 6); // Comb. PID / Bachelor+Baryon+Meson
2818         }
2819
2820         lContainerPIDVars[0] = 0.; lContainerPIDVars[1] = 0.; lContainerPIDVars[2] = 0.; lContainerPIDVars[3] = 0.;
2821
2822         // Xi Plus
2823         if( lChargeXi > 0 && lAssoXiPlus ) {
2824                 lContainerPIDVars[0] = lmcPt              ;
2825                 lContainerPIDVars[1] = lInvMassXiPlus     ;
2826                 lContainerPIDVars[2] = lmcRapCasc           ;
2827                 if (fCollidingSystems) lContainerPIDVars[3] = lcentrality;
2828                 else lContainerPIDVars[3] = lSPDTrackletsMultiplicity ;
2829
2830                 // No PID
2831                         fCFContCascadePIDAsXiPlus->Fill(lContainerPIDVars, 0); // No PID
2832                 // TPC PID
2833                 if( lIsBachelorPionForTPC  )
2834                         fCFContCascadePIDAsXiPlus->Fill(lContainerPIDVars, 1); // TPC PID / 4-#sigma cut on Bachelor track
2835
2836                 if( lIsBachelorPionForTPC &&
2837                     lIsNegProtonForTPC     )
2838                         fCFContCascadePIDAsXiPlus->Fill(lContainerPIDVars, 2); // TPC PID / 4-#sigma cut on Bachelor+Baryon tracks
2839
2840                 if( lIsBachelorPionForTPC &&
2841                     lIsNegProtonForTPC    &&
2842                     lIsPosPionForTPC       )
2843                         fCFContCascadePIDAsXiPlus->Fill(lContainerPIDVars, 3); // TPC PID / 4-#sigma cut on Bachelor+Baryon+Meson tracks
2844
2845                 // Combined PID
2846                 if( lIsBachelorPion        )
2847                         fCFContCascadePIDAsXiPlus->Fill(lContainerPIDVars, 4); // Comb. PID / Bachelor
2848
2849                 if( lIsBachelorPion       &&
2850                     lIsNegInXiProton    )
2851                         fCFContCascadePIDAsXiPlus->Fill(lContainerPIDVars, 5); // Comb. PID / Bachelor+Baryon
2852
2853                 if(lIsBachelorPion     &&
2854                    lIsNegInXiProton &&
2855                    lIsPosInXiPion    )
2856                         fCFContCascadePIDAsXiPlus->Fill(lContainerPIDVars, 6); // Comb. PID / Bachelor+Baryon+Meson
2857         }
2858
2859         lContainerPIDVars[0] = 0.; lContainerPIDVars[1] = 0.; lContainerPIDVars[2] = 0.; lContainerPIDVars[3] = 0.;
2860
2861         // Omega Minus
2862         if( lChargeXi < 0 && lAssoOmegaMinus ) {
2863                 lContainerPIDVars[0] = lmcPt              ;
2864                 lContainerPIDVars[1] = lInvMassOmegaMinus ;
2865                 lContainerPIDVars[2] = lmcRapCasc         ;
2866                 if (fCollidingSystems) lContainerPIDVars[3] = lcentrality;
2867                 else lContainerPIDVars[3] = lSPDTrackletsMultiplicity ;
2868
2869                 // No PID
2870                         fCFContCascadePIDAsOmegaMinus->Fill(lContainerPIDVars, 0); // No PID
2871                 // TPC PID
2872                 if( lIsBachelorKaonForTPC  )
2873                         fCFContCascadePIDAsOmegaMinus->Fill(lContainerPIDVars, 1); // TPC PID / 4-#sigma cut on Bachelor track
2874
2875                 if( lIsBachelorKaonForTPC &&
2876                     lIsPosProtonForTPC     )
2877                         fCFContCascadePIDAsOmegaMinus->Fill(lContainerPIDVars, 2); // TPC PID / 4-#sigma cut on Bachelor+Baryon tracks
2878
2879                 if( lIsBachelorKaonForTPC &&
2880                     lIsPosProtonForTPC    &&
2881                     lIsNegPionForTPC       )
2882                         fCFContCascadePIDAsOmegaMinus->Fill(lContainerPIDVars, 3); // TPC PID / 4-#sigma cut on Bachelor+Baryon+Meson tracks
2883
2884                 // Combined PID
2885                 if( lIsBachelorKaon        )
2886                         fCFContCascadePIDAsOmegaMinus->Fill(lContainerPIDVars, 4); // Comb. PID / Bachelor
2887
2888                 if( lIsBachelorKaon       &&
2889                     lIsPosInOmegaProton    )
2890                         fCFContCascadePIDAsOmegaMinus->Fill(lContainerPIDVars, 5); // Comb. PID / Bachelor+Baryon
2891
2892                 if(lIsBachelorKaon     &&
2893                    lIsPosInOmegaProton &&
2894                    lIsNegInOmegaPion    )
2895                         fCFContCascadePIDAsOmegaMinus->Fill(lContainerPIDVars, 6); // Comb. PID / Bachelor+Baryon+Meson
2896         }
2897
2898         lContainerPIDVars[0] = 0.; lContainerPIDVars[1] = 0.; lContainerPIDVars[2] = 0.; lContainerPIDVars[3] = 0.;
2899
2900         // Omega Plus
2901         if( lChargeXi > 0 && lAssoOmegaPlus) {
2902                 lContainerPIDVars[0] = lmcPt              ;
2903                 lContainerPIDVars[1] = lInvMassOmegaPlus  ;
2904                 lContainerPIDVars[2] = lmcRapCasc         ;
2905                 if (fCollidingSystems) lContainerPIDVars[3] = lcentrality;
2906                 else lContainerPIDVars[3] = lSPDTrackletsMultiplicity ;
2907
2908                 // No PID
2909                         fCFContCascadePIDAsOmegaPlus->Fill(lContainerPIDVars, 0); // No PID
2910                 // TPC PID
2911                 if( lIsBachelorKaonForTPC  )
2912                         fCFContCascadePIDAsOmegaPlus->Fill(lContainerPIDVars, 1); // TPC PID / 4-#sigma cut on Bachelor track
2913
2914                 if( lIsBachelorKaonForTPC &&
2915                     lIsNegProtonForTPC     )
2916                         fCFContCascadePIDAsOmegaPlus->Fill(lContainerPIDVars, 2); // TPC PID / 4-#sigma cut on Bachelor+Baryon tracks
2917
2918                 if( lIsBachelorKaonForTPC &&
2919                     lIsNegProtonForTPC    &&
2920                     lIsPosPionForTPC       )
2921                         fCFContCascadePIDAsOmegaPlus->Fill(lContainerPIDVars, 3); // TPC PID / 4-#sigma cut on Bachelor+Baryon+Meson tracks
2922
2923                 // Combined PID
2924                 if( lIsBachelorKaon        )
2925                         fCFContCascadePIDAsOmegaPlus->Fill(lContainerPIDVars, 4); // Comb. PID / Bachelor
2926
2927                 if( lIsBachelorKaon       &&
2928                     lIsNegInOmegaProton    )
2929                         fCFContCascadePIDAsOmegaPlus->Fill(lContainerPIDVars, 5); // Comb. PID / Bachelor+Baryon
2930
2931                 if(lIsBachelorKaon     &&
2932                    lIsNegInOmegaProton &&
2933                    lIsPosInOmegaPion    )
2934                         fCFContCascadePIDAsOmegaPlus->Fill(lContainerPIDVars, 6); // Comb. PID / Bachelor+Baryon+Meson
2935         }
2936
2937       }
2938
2939 }// End of loop over reconstructed cascades
2940
2941 fHistnAssoXiMinus->Fill(nAssoXiMinus);
2942 fHistnAssoXiPlus->Fill(nAssoXiPlus);
2943 fHistnAssoOmegaMinus->Fill(nAssoOmegaMinus);
2944 fHistnAssoOmegaPlus->Fill(nAssoOmegaPlus);
2945 //Printf("N asso Xi- = %n ", nAssoXiMinus);
2946
2947   // Post output data.
2948  PostData(1, fListHistCascade);
2949 }
2950
2951
2952 Int_t AliAnalysisTaskCheckPerformanceCascadePbPb::DoESDTrackWithTPCrefitMultiplicity(const AliESDEvent *lESDevent) {
2953     // Checking the number of tracks with TPCrefit for each event
2954     // Needed for a rough assessment of the event multiplicity
2955
2956         Int_t nTrackWithTPCrefitMultiplicity = 0;
2957         for (Int_t iTrackIdx = 0; iTrackIdx < (InputEvent())->GetNumberOfTracks(); iTrackIdx++) {
2958                 AliESDtrack *esdTrack   = 0x0;
2959                              esdTrack   = lESDevent->GetTrack( iTrackIdx );
2960                 if (!esdTrack) { AliWarning("Pb / Could not retrieve one track within the track loop for TPCrefit check ..."); continue; }
2961
2962                 ULong_t lTrackStatus    = esdTrack->GetStatus();
2963                             if ((lTrackStatus&AliESDtrack::kTPCrefit)    == 0) continue;
2964                             else nTrackWithTPCrefitMultiplicity++;
2965                     // FIXME :
2966                     // The goal here is to get a better assessment of the event multiplicity.
2967                     // (InputEvent())->GetNumberOfTracks() takes into account ITS std alone tracks + global tracks
2968                     // This may introduce a bias. Hence the number of TPC refit tracks.
2969                     // Note : the event multiplicity = analysis on its own... See Jacek's or Jan Fiete's analysis on dN/d(pt) and dN/d(eta)
2970
2971         }// end loop over all event tracks
2972         return  nTrackWithTPCrefitMultiplicity;
2973 }
2974
2975 //________________________________________________________________________
2976 void AliAnalysisTaskCheckPerformanceCascadePbPb::Terminate(Option_t *) {
2977   // Draw result to the screen
2978   // Called once at the end of the query
2979
2980   TList *cRetrievedList = 0x0;
2981   cRetrievedList = (TList*)GetOutputData(1);
2982   if(!cRetrievedList) {
2983         Printf("ERROR - AliAnalysisTaskCheckPerformanceCascadePbPb : ouput data container list not available\n");
2984         return;
2985   }
2986
2987   fHistMCTrackMultiplicity = dynamic_cast<TH1F*> (  cRetrievedList->FindObject("fHistMCTrackMultiplicity")  );
2988   if (!fHistMCTrackMultiplicity) {
2989     Printf("ERROR - AliAnalysisTaskCheckPerformanceCascadePbPb : fHistMCTrackMultiplicity not available");
2990     return;
2991   }
2992
2993
2994   TCanvas *canCheckPerformanceCascade = new TCanvas("AliAnalysisTaskCheckPerformanceCascadePbPb","Multiplicity",10,10,510,510);
2995   canCheckPerformanceCascade->cd(1)->SetLogy();
2996
2997   fHistMCTrackMultiplicity->SetMarkerStyle(22);
2998   fHistMCTrackMultiplicity->DrawCopy("E");
2999
3000 }