PWGLF/STRANGENESS/Cascades/AliAnalysisTaskCheckPerformanceCascadepp276.cxx

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