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