PWG2/SPECTRA/AliAnalysisTaskCheckCascade.cxx

   1
   2 /**************************************************************************
   3  *  Authors : Antonin Maire, Boris Hippolyte                              *
   4  * Contributors are mentioned in the code where appropriate.              *
   5  *                                                                        *
   6  * Permission to use, copy, modify and distribute this software and its   *
   7  * documentation strictly for non-commercial purposes is hereby granted   *
   8  * without fee, provided that the above copyright notice appears in all   *
   9  * copies and that both the copyright notice and this permission notice   *
  10  * appear in the supporting documentation. The authors make no claims     *
  11  * about the suitability of this software for any purpose. It is          *
  12  * provided "as is" without express or implied warranty.                  *
  13  **************************************************************************/
  14
  15 //-----------------------------------------------------------------
  16 //                 AliAnalysisTaskCheckCascade class
  17 //            (AliAnalysisTaskCheckCascade)
  18 //            This task has four roles :
  19 //              1. QAing the Cascades from ESD and AOD
  20 //                 Origin:  AliAnalysisTaskESDCheckV0 by B.H. Nov2007, hippolyt@in2p3.fr
  21 //              2. Prepare the plots which stand as raw material for yield extraction (wi/wo PID)
  22 //              3. Supply an AliCFContainer meant to define the optimised topological selections
  23 //              4. Rough azimuthal correlation study (Eta, Phi)
  24 //            Adapted to Cascade : A.Maire Mar2008, antonin.maire@ires.in2p3.fr
  25 //            Modified :           A.Maire Nov2009, antonin.maire@ires.in2p3.fr
  26 //-----------------------------------------------------------------
  27
  28
  29
  30 class TTree;
  31 class TParticle;
  32 class TVector3;
  33
  34 //class AliMCEventHandler;
  35 //class AliMCEvent;
  36 //class AliStack;
  37
  38 class AliESDVertex;
  39 class AliAODVertex;
  40 class AliESDv0;
  41 class AliAODv0;
  42
  43 #include <Riostream.h>
  44
  45 #include "TList.h"
  46 #include "TH1.h"
  47 #include "TH2.h"
  48 #include "TH3.h"
  49 #include "THnSparse.h"
  50 #include "TVector3.h"
  51 #include "TCanvas.h"
  52 #include "TMath.h"
  53
  54
  55 #include "AliLog.h"
  56
  57 #include "AliESDEvent.h"
  58 #include "AliAODEvent.h"
  59 //#include "AliCascadeVertexer.h"
  60 #include "AliTPCpidESD.h"
  61 #include "AliCFContainer.h"
  62 #include "AliMultiplicity.h"
  63
  64 #include "AliESDcascade.h"
  65 #include "AliAODcascade.h"
  66
  67 #include "AliAnalysisTaskCheckCascade.h"
  68
  69 ClassImp(AliAnalysisTaskCheckCascade)
  70
  71
  72
  73 //________________________________________________________________________
  74 AliAnalysisTaskCheckCascade::AliAnalysisTaskCheckCascade()
  75   : AliAnalysisTaskSE(), fAnalysisType("ESD"), fCollidingSystems(0), fTpcPidManager(0),
  76
  77         // - Cascade part initialisation
  78     fListHistCascade(0),
  79     fHistTrackMultiplicity(0), fHistCascadeMultiplicity(0),
  80     fHistVtxStatus(0),
  81
  82     fHistPosTrkgPrimaryVtxX(0), fHistPosTrkgPrimaryVtxY(0), fHistPosTrkgPrimaryVtxZ(0), fHistTrkgPrimaryVtxRadius(0),
  83     fHistPosBestPrimaryVtxX(0), fHistPosBestPrimaryVtxY(0), fHistPosBestPrimaryVtxZ(0), fHistBestPrimaryVtxRadius(0),
  84     f2dHistTrkgPrimVtxVsBestPrimVtx(0),
  85
  86     fHistEffMassXi(0),  fHistChi2Xi(0),
  87     fHistDcaXiDaughters(0), fHistDcaBachToPrimVertex(0), fHistXiCosineOfPointingAngle(0), fHistXiRadius(0),
  88
  89     fHistMassLambdaAsCascDghter(0),
  90     fHistV0Chi2Xi(0),
  91     fHistDcaV0DaughtersXi(0),
  92     fHistDcaV0ToPrimVertexXi(0),
  93     fHistV0CosineOfPointingAngleXi(0),
  94     fHistV0RadiusXi(0),
  95     fHistDcaPosToPrimVertexXi(0), fHistDcaNegToPrimVertexXi(0),
  96
  97     fHistMassXiMinus(0), fHistMassXiPlus(0),
  98     fHistMassOmegaMinus(0), fHistMassOmegaPlus(0),
  99     fHistMassWithCombPIDXiMinus(0), fHistMassWithCombPIDXiPlus(0),
 100     fHistMassWithCombPIDOmegaMinus(0), fHistMassWithCombPIDOmegaPlus(0),
 101
 102     fHistXiTransvMom(0),    fHistXiTotMom(0),
 103     fHistBachTransvMom(0),   fHistBachTotMom(0),
 104
 105     fHistChargeXi(0),
 106     fHistV0toXiCosineOfPointingAngle(0),
 107
 108     fHistRapXi(0), fHistRapOmega(0), fHistEta(0),
 109     fHistTheta(0), fHistPhi(0),
 110
 111     f2dHistArmenteros(0),
 112     f2dHistEffMassLambdaVsEffMassXiMinus(0), f2dHistEffMassXiVsEffMassOmegaMinus(0),
 113     f2dHistEffMassLambdaVsEffMassXiPlus(0), f2dHistEffMassXiVsEffMassOmegaPlus(0),
 114     f2dHistXiRadiusVsEffMassXiMinus(0), f2dHistXiRadiusVsEffMassXiPlus(0),
 115     f2dHistXiRadiusVsEffMassOmegaMinus(0), f2dHistXiRadiusVsEffMassOmegaPlus(0),
 116
 117     f3dHistXiPtVsEffMassVsYXiMinus(0), f3dHistXiPtVsEffMassVsYXiPlus(0),
 118     f3dHistXiPtVsEffMassVsYOmegaMinus(0), f3dHistXiPtVsEffMassVsYOmegaPlus(0),
 119
 120     f3dHistXiPtVsEffMassVsYWithCombPIDXiMinus(0), f3dHistXiPtVsEffMassVsYWithCombPIDXiPlus(0),
 121     f3dHistXiPtVsEffMassVsYWithCombPIDOmegaMinus(0),  f3dHistXiPtVsEffMassVsYWithCombPIDOmegaPlus(0),
 122     f3dHistXiPtVsEffMassVsYWith2CombPIDOmegaMinus(0), f3dHistXiPtVsEffMassVsYWith2CombPIDOmegaPlus(0),
 123     f3dHistXiPtVsEffMassVsYWithTpcPIDOmegaMinus(0),
 124
 125     fCFContCascadePIDXiMinus(0),
 126     fCFContCascadePIDXiPlus(0),
 127     fCFContCascadePIDOmegaMinus(0),
 128     fCFContCascadePIDOmegaPlus(0),
 129     fCFContCascadeCuts(0),
 130
 131     fHnSpAngularCorrXiMinus(0), fHnSpAngularCorrXiPlus(0),
 132     fHnSpAngularCorrOmegaMinus(0), fHnSpAngularCorrOmegaPlus(0)
 133
 134 {
 135   // Dummy Constructor
 136 }
 137
 138
 139
 140
 141
 142
 143
 144
 145 //________________________________________________________________________
 146 AliAnalysisTaskCheckCascade::AliAnalysisTaskCheckCascade(const char *name)
 147   : AliAnalysisTaskSE(name), fAnalysisType("ESD"), fCollidingSystems(0), fTpcPidManager(0),
 148
 149         // - Cascade part initialisation
 150     fListHistCascade(0),
 151     fHistTrackMultiplicity(0), fHistCascadeMultiplicity(0),
 152     fHistVtxStatus(0),
 153
 154     fHistPosTrkgPrimaryVtxX(0), fHistPosTrkgPrimaryVtxY(0), fHistPosTrkgPrimaryVtxZ(0), fHistTrkgPrimaryVtxRadius(0),
 155     fHistPosBestPrimaryVtxX(0), fHistPosBestPrimaryVtxY(0), fHistPosBestPrimaryVtxZ(0), fHistBestPrimaryVtxRadius(0),
 156     f2dHistTrkgPrimVtxVsBestPrimVtx(0),
 157
 158     fHistEffMassXi(0),  fHistChi2Xi(0),
 159     fHistDcaXiDaughters(0), fHistDcaBachToPrimVertex(0), fHistXiCosineOfPointingAngle(0), fHistXiRadius(0),
 160
 161     fHistMassLambdaAsCascDghter(0),
 162     fHistV0Chi2Xi(0),
 163     fHistDcaV0DaughtersXi(0),
 164     fHistDcaV0ToPrimVertexXi(0),
 165     fHistV0CosineOfPointingAngleXi(0),
 166     fHistV0RadiusXi(0),
 167     fHistDcaPosToPrimVertexXi(0), fHistDcaNegToPrimVertexXi(0),
 168
 169     fHistMassXiMinus(0), fHistMassXiPlus(0),
 170     fHistMassOmegaMinus(0), fHistMassOmegaPlus(0),
 171     fHistMassWithCombPIDXiMinus(0), fHistMassWithCombPIDXiPlus(0),
 172     fHistMassWithCombPIDOmegaMinus(0), fHistMassWithCombPIDOmegaPlus(0),
 173
 174     fHistXiTransvMom(0),    fHistXiTotMom(0),
 175     fHistBachTransvMom(0),   fHistBachTotMom(0),
 176
 177     fHistChargeXi(0),
 178     fHistV0toXiCosineOfPointingAngle(0),
 179
 180     fHistRapXi(0), fHistRapOmega(0), fHistEta(0),
 181     fHistTheta(0), fHistPhi(0),
 182
 183     f2dHistArmenteros(0),
 184     f2dHistEffMassLambdaVsEffMassXiMinus(0), f2dHistEffMassXiVsEffMassOmegaMinus(0),
 185     f2dHistEffMassLambdaVsEffMassXiPlus(0), f2dHistEffMassXiVsEffMassOmegaPlus(0),
 186     f2dHistXiRadiusVsEffMassXiMinus(0), f2dHistXiRadiusVsEffMassXiPlus(0),
 187     f2dHistXiRadiusVsEffMassOmegaMinus(0), f2dHistXiRadiusVsEffMassOmegaPlus(0),
 188
 189     f3dHistXiPtVsEffMassVsYXiMinus(0), f3dHistXiPtVsEffMassVsYXiPlus(0),
 190     f3dHistXiPtVsEffMassVsYOmegaMinus(0), f3dHistXiPtVsEffMassVsYOmegaPlus(0),
 191
 192     f3dHistXiPtVsEffMassVsYWithCombPIDXiMinus(0), f3dHistXiPtVsEffMassVsYWithCombPIDXiPlus(0),
 193     f3dHistXiPtVsEffMassVsYWithCombPIDOmegaMinus(0),  f3dHistXiPtVsEffMassVsYWithCombPIDOmegaPlus(0),
 194     f3dHistXiPtVsEffMassVsYWith2CombPIDOmegaMinus(0), f3dHistXiPtVsEffMassVsYWith2CombPIDOmegaPlus(0),
 195     f3dHistXiPtVsEffMassVsYWithTpcPIDOmegaMinus(0),
 196
 197     fCFContCascadePIDXiMinus(0),
 198     fCFContCascadePIDXiPlus(0),
 199     fCFContCascadePIDOmegaMinus(0),
 200     fCFContCascadePIDOmegaPlus(0),
 201     fCFContCascadeCuts(0),
 202
 203     fHnSpAngularCorrXiMinus(0), fHnSpAngularCorrXiPlus(0),
 204     fHnSpAngularCorrOmegaMinus(0), fHnSpAngularCorrOmegaPlus(0)
 205
 206 {
 207   // Constructor
 208
 209   // Define input and output slots here
 210   // Input slot #0 works with a TChain
 211
 212   // Output slot #0 writes into a TList container (Cascade)
 213   DefineOutput(1, TList::Class());
 214 }
 215
 216
 217
 218
 219
 220
 221 //________________________________________________________________________
 222 void AliAnalysisTaskCheckCascade::UserCreateOutputObjects()
 223 {
 224   // Create histograms
 225   // Called once
 226
 227
 228
 229  fListHistCascade = new TList();
 230
 231
 232         // - General histos
 233
 234 if(! fHistTrackMultiplicity) {
 235         if(fCollidingSystems)// AA collisions
 236                 fHistTrackMultiplicity = new TH1F("fHistTrackMultiplicity",
 237                         "Multiplicity distribution;Number of tracks;Events",
 238                         200, 0, 20000);
 239         else // pp collisions
 240                 fHistTrackMultiplicity = new TH1F("fHistTrackMultiplicity",
 241                         "Track Multiplicity;Nbr of tracks/Evt;Events",
 242                         200, 0, 200);
 243         fListHistCascade->Add(fHistTrackMultiplicity);
 244 }
 245
 246 if(! fHistCascadeMultiplicity) {
 247         if(fCollidingSystems)// AA collisions
 248                 fHistCascadeMultiplicity = new TH1F("fHistCascadeMultiplicity",
 249                         "Multiplicity distribution;Number of Cascades;Events",
 250                         50, 0, 50);
 251         else // pp collisions
 252                 fHistCascadeMultiplicity = new TH1F("fHistCascadeMultiplicity",
 253                         "Cascades per event;Nbr of Cascades/Evt;Events",
 254                         10, 0, 10);
 255         fListHistCascade->Add(fHistCascadeMultiplicity);
 256 }
 257
 258
 259
 260 if(! fHistVtxStatus ){
 261         fHistVtxStatus   = new TH1F( "fHistVtxStatus" , "Does a Trckg Prim.vtx exist ?; true=1 or false=0; Nb of Events" , 4, -1.0, 3.0 );
 262         fListHistCascade->Add(fHistVtxStatus);
 263 }
 264
 265
 266
 267
 268
 269
 270         // - Vertex Positions
 271
 272 if(! fHistPosTrkgPrimaryVtxX ){
 273         fHistPosTrkgPrimaryVtxX   = new TH1F( "fHistPosTrkgPrimaryVtxX" , "Trkg Prim. Vertex Position in x; x (cm); Events" , 200, -0.5, 0.5 );
 274         fListHistCascade->Add(fHistPosTrkgPrimaryVtxX);
 275 }
 276
 277
 278 if(! fHistPosTrkgPrimaryVtxY){
 279         fHistPosTrkgPrimaryVtxY   = new TH1F( "fHistPosTrkgPrimaryVtxY" , "Trkg Prim. Vertex Position in y; y (cm); Events" , 200, -0.5, 0.5 );
 280         fListHistCascade->Add(fHistPosTrkgPrimaryVtxY);
 281 }
 282
 283 if(! fHistPosTrkgPrimaryVtxZ ){
 284         fHistPosTrkgPrimaryVtxZ   = new TH1F( "fHistPosTrkgPrimaryVtxZ" , "Trkg Prim. Vertex Position in z; z (cm); Events" , 100, -15.0, 15.0 );
 285         fListHistCascade->Add(fHistPosTrkgPrimaryVtxZ);
 286 }
 287
 288 if(! fHistTrkgPrimaryVtxRadius ){
 289         fHistTrkgPrimaryVtxRadius  = new TH1F( "fHistTrkgPrimaryVtxRadius",  "Trkg Prim. Vertex radius; r (cm); Events" , 150, 0., 15.0 );
 290         fListHistCascade->Add(fHistTrkgPrimaryVtxRadius);
 291 }
 292
 293
 294
 295
 296 if(! fHistPosBestPrimaryVtxX ){
 297         fHistPosBestPrimaryVtxX   = new TH1F( "fHistPosBestPrimaryVtxX" , "Best Prim. Vertex Position in x; x (cm); Events" , 200, -0.5, 0.5 );
 298         fListHistCascade->Add(fHistPosBestPrimaryVtxX);
 299 }
 300
 301 if(! fHistPosBestPrimaryVtxY){
 302         fHistPosBestPrimaryVtxY   = new TH1F( "fHistPosBestPrimaryVtxY" , "Best Prim. Vertex Position in y; y (cm); Events" , 200, -0.5, 0.5 );
 303         fListHistCascade->Add(fHistPosBestPrimaryVtxY);
 304 }
 305
 306 if(! fHistPosBestPrimaryVtxZ ){
 307         fHistPosBestPrimaryVtxZ   = new TH1F( "fHistPosBestPrimaryVtxZ" , "Best Prim. Vertex Position in z; z (cm); Events" , 100, -15.0, 15.0 );
 308         fListHistCascade->Add(fHistPosBestPrimaryVtxZ);
 309 }
 310
 311 if(! fHistBestPrimaryVtxRadius ){
 312         fHistBestPrimaryVtxRadius  = new TH1F( "fHistBestPrimaryVtxRadius",  "Best Prim.  vertex radius; r (cm); Events" , 150, 0., 15.0 );
 313         fListHistCascade->Add(fHistBestPrimaryVtxRadius);
 314 }
 315
 316 if(! f2dHistTrkgPrimVtxVsBestPrimVtx) {
 317         f2dHistTrkgPrimVtxVsBestPrimVtx = new TH2F( "f2dHistTrkgPrimVtxVsBestPrimVtx", "r_{Trck Prim. Vtx} Vs r_{Best Prim. Vtx}; r_{Track Vtx} (cm); r_{Best Vtx} (cm)", 300, 0., 15.0, 300, 0., 15.);
 318         fListHistCascade->Add(f2dHistTrkgPrimVtxVsBestPrimVtx);
 319 }
 320
 321
 322
 323
 324 // - Typical histos for cascades
 325
 326
 327 if(! fHistEffMassXi) {
 328      fHistEffMassXi = new TH1F("fHistEffMassXi", "Cascade candidates ; Invariant Mass (GeV/c^{2}) ; Counts", 200, 1.2, 2.0);
 329      fListHistCascade->Add(fHistEffMassXi);
 330 }
 331
 332 if(! fHistChi2Xi ){
 333         fHistChi2Xi = new TH1F("fHistChi2Xi", "Cascade #chi^{2}; #chi^{2}; Number of Cascades", 160, 0, 40);
 334         fListHistCascade->Add(fHistChi2Xi);
 335 }
 336
 337 if(! fHistDcaXiDaughters ){
 338         fHistDcaXiDaughters = new TH1F( "fHistDcaXiDaughters",  "DCA between Xi Daughters; DCA (cm) ; Number of Cascades", 100, 0., 0.5);
 339         fListHistCascade->Add(fHistDcaXiDaughters);
 340 }
 341
 342 if(! fHistDcaBachToPrimVertex) {
 343         fHistDcaBachToPrimVertex = new TH1F("fHistDcaBachToPrimVertex", "DCA of Bach. to Prim. Vertex;DCA (cm);Number of Cascades", 250, 0., 0.25);
 344         fListHistCascade->Add(fHistDcaBachToPrimVertex);
 345 }
 346
 347 if(! fHistXiCosineOfPointingAngle) {
 348         fHistXiCosineOfPointingAngle = new TH1F("fHistXiCosineOfPointingAngle", "Cosine of Xi Pointing Angle; Cos (Xi Point.Angl);Number of Xis", 200, 0.99, 1.0);
 349         fListHistCascade->Add(fHistXiCosineOfPointingAngle);
 350 }
 351
 352 if(! fHistXiRadius ){
 353         fHistXiRadius  = new TH1F( "fHistXiRadius",  "Casc. decay transv. radius; r (cm); Counts" , 1050, 0., 105.0 );
 354         fListHistCascade->Add(fHistXiRadius);
 355 }
 356
 357
 358 // - Histos about ~ the "V0 part" of the cascade,  coming by inheritance from AliESDv0
 359
 360
 361
 362 if (! fHistMassLambdaAsCascDghter) {
 363      fHistMassLambdaAsCascDghter = new TH1F("fHistMassLambdaAsCascDghter","#Lambda associated to Casc. candidates;Eff. Mass (GeV/c^{2});Counts", 300,1.00,1.3);
 364     fListHistCascade->Add(fHistMassLambdaAsCascDghter);
 365 }
 366
 367 if (! fHistV0Chi2Xi) {
 368         fHistV0Chi2Xi = new TH1F("fHistV0Chi2Xi", "V0 #chi^{2}, in cascade; #chi^{2};Counts", 160, 0, 40);
 369         fListHistCascade->Add(fHistV0Chi2Xi);
 370 }
 371
 372 if (! fHistDcaV0DaughtersXi) {
 373         fHistDcaV0DaughtersXi = new TH1F("fHistDcaV0DaughtersXi", "DCA between V0 daughters, in cascade;DCA (cm);Number of V0s", 120, 0., 0.6);
 374         fListHistCascade->Add(fHistDcaV0DaughtersXi);
 375 }
 376
 377 if (! fHistDcaV0ToPrimVertexXi) {
 378         fHistDcaV0ToPrimVertexXi = new TH1F("fHistDcaV0ToPrimVertexXi", "DCA of V0 to Prim. Vertex, in cascade;DCA (cm);Number of Cascades", 200, 0., 1.);
 379         fListHistCascade->Add(fHistDcaV0ToPrimVertexXi);
 380 }
 381
 382 if (! fHistV0CosineOfPointingAngleXi) {
 383         fHistV0CosineOfPointingAngleXi = new TH1F("fHistV0CosineOfPointingAngleXi", "Cosine of V0 Pointing Angle, in cascade;Cos(V0 Point. Angl); Counts", 200, 0.98, 1.0);
 384         fListHistCascade->Add(fHistV0CosineOfPointingAngleXi);
 385 }
 386
 387 if (! fHistV0RadiusXi) {
 388         fHistV0RadiusXi = new TH1F("fHistV0RadiusXi", "V0 decay radius, in cascade; radius (cm); Counts", 1050, 0., 105.0);
 389         fListHistCascade->Add(fHistV0RadiusXi);
 390 }
 391
 392 if (! fHistDcaPosToPrimVertexXi) {
 393         fHistDcaPosToPrimVertexXi = new TH1F("fHistDcaPosToPrimVertexXi", "DCA of V0 pos daughter to Prim. Vertex;DCA (cm);Counts", 300, 0, 3);
 394         fListHistCascade->Add(fHistDcaPosToPrimVertexXi);
 395 }
 396
 397 if (! fHistDcaNegToPrimVertexXi) {
 398         fHistDcaNegToPrimVertexXi = new TH1F("fHistDcaNegToPrimVertexXi", "DCA of V0 neg daughter to Prim. Vertex;DCA (cm);Counts", 300, 0, 3);
 399         fListHistCascade->Add(fHistDcaNegToPrimVertexXi);
 400 }
 401
 402
 403
 404
 405         // - Effective mass histos for cascades.
 406 // By cascade hyp
 407 if (! fHistMassXiMinus) {
 408     fHistMassXiMinus = new TH1F("fHistMassXiMinus","#Xi^{-} candidates;M( #Lambda , #pi^{-} ) (GeV/c^{2});Counts", 200,1.2,2.0);
 409     fListHistCascade->Add(fHistMassXiMinus);
 410 }
 411
 412 if (! fHistMassXiPlus) {
 413     fHistMassXiPlus = new TH1F("fHistMassXiPlus","#Xi^{+} candidates;M( #bar{#Lambda}^{0} , #pi^{+} ) (GeV/c^{2});Counts",200,1.2,2.0);
 414     fListHistCascade->Add(fHistMassXiPlus);
 415 }
 416
 417 if (! fHistMassOmegaMinus) {
 418         fHistMassOmegaMinus = new TH1F("fHistMassOmegaMinus","#Omega^{-} candidates;M( #Lambda , K^{-} ) (GeV/c^{2});Counts", 250,1.5,2.5);
 419     fListHistCascade->Add(fHistMassOmegaMinus);
 420 }
 421
 422 if (! fHistMassOmegaPlus) {
 423         fHistMassOmegaPlus = new TH1F("fHistMassOmegaPlus","#Omega^{+} candidates;M( #bar{#Lambda}^{0} , K^{+} ) (GeV/c^{2});Counts", 250,1.5,2.5);
 424     fListHistCascade->Add(fHistMassOmegaPlus);
 425 }
 426
 427 // By cascade hyp + bachelor PID
 428 if (! fHistMassWithCombPIDXiMinus) {
 429     fHistMassWithCombPIDXiMinus = new TH1F("fHistMassWithCombPIDXiMinus","#Xi^{-} candidates, with Bach. comb. PID;M( #Lambda , #pi^{-} ) (GeV/c^{2});Counts", 200,1.2,2.0);
 430     fListHistCascade->Add(fHistMassWithCombPIDXiMinus);
 431 }
 432
 433 if (! fHistMassWithCombPIDXiPlus) {
 434     fHistMassWithCombPIDXiPlus = new TH1F("fHistMassWithCombPIDXiPlus","#Xi^{+} candidates, with Bach. comb. PID;M( #bar{#Lambda}^{0} , #pi^{+} ) (GeV/c^{2});Counts",200,1.2,2.0);
 435     fListHistCascade->Add(fHistMassWithCombPIDXiPlus);
 436 }
 437
 438 if (! fHistMassWithCombPIDOmegaMinus) {
 439         fHistMassWithCombPIDOmegaMinus = new TH1F("fHistMassWithCombPIDOmegaMinus","#Omega^{-} candidates, with Bach. comb. PID;M( #Lambda , K^{-} ) (GeV/c^{2});Counts", 250,1.5,2.5);
 440     fListHistCascade->Add(fHistMassWithCombPIDOmegaMinus);
 441 }
 442
 443 if (! fHistMassWithCombPIDOmegaPlus) {
 444         fHistMassWithCombPIDOmegaPlus = new TH1F("fHistMassWithCombPIDOmegaPlus","#Omega^{+} candidates, with Bach. comb. PID;M( #bar{#Lambda}^{0} , K^{+} ) (GeV/c^{2});Counts", 250,1.5,2.5);
 445     fListHistCascade->Add(fHistMassWithCombPIDOmegaPlus);
 446 }
 447
 448
 449
 450         // - Complements for QA
 451
 452 if(! fHistXiTransvMom ){
 453         fHistXiTransvMom  = new TH1F( "fHistXiTransvMom" , "Xi transverse momentum ; p_{t}(#Xi) (GeV/c); Counts", 100, 0.0, 10.0);
 454         fListHistCascade->Add(fHistXiTransvMom);
 455 }
 456
 457 if(! fHistXiTotMom ){
 458         fHistXiTotMom  = new TH1F( "fHistXiTotMom" , "Xi momentum norm; p_{tot}(#Xi) (GeV/c); Counts", 150, 0.0, 15.0);
 459         fListHistCascade->Add(fHistXiTotMom);
 460 }
 461
 462
 463 if(! fHistBachTransvMom ){
 464         fHistBachTransvMom  = new TH1F( "fHistBachTransvMom" , "Bach. transverse momentum ; p_{t}(Bach.) (GeV/c); Counts", 100, 0.0, 5.0);
 465         fListHistCascade->Add(fHistBachTransvMom);
 466 }
 467
 468 if(! fHistBachTotMom ){
 469         fHistBachTotMom  = new TH1F( "fHistBachTotMom" , "Bach. momentum norm; p_{tot}(Bach.) (GeV/c); Counts", 100, 0.0, 5.0);
 470         fListHistCascade->Add(fHistBachTotMom);
 471 }
 472
 473
 474 if(! fHistChargeXi ){
 475         fHistChargeXi  = new TH1F( "fHistChargeXi" , "Charge of casc. candidates ; Sign ; Counts", 5, -2.0, 3.0);
 476         fListHistCascade->Add(fHistChargeXi);
 477 }
 478
 479
 480 if (! fHistV0toXiCosineOfPointingAngle) {
 481         fHistV0toXiCosineOfPointingAngle = new TH1F("fHistV0toXiCosineOfPointingAngle", "Cos. of V0 Ptng Angl / Xi vtx ;Cos(V0 Point. Angl / Xi vtx); Counts", 100, 0.99, 1.0);
 482         fListHistCascade->Add(fHistV0toXiCosineOfPointingAngle);
 483 }
 484
 485
 486 if(! fHistRapXi ){
 487         fHistRapXi  = new TH1F( "fHistRapXi" , "Rapidity of Xi candidates ; y ; Counts", 200, -5.0, 5.0);
 488         fListHistCascade->Add(fHistRapXi);
 489 }
 490
 491 if(! fHistRapOmega ){
 492         fHistRapOmega  = new TH1F( "fHistRapOmega" , "Rapidity of Omega candidates ; y ; Counts", 200, -5.0, 5.0);
 493         fListHistCascade->Add(fHistRapOmega);
 494 }
 495
 496 if(! fHistEta ){
 497         fHistEta  = new TH1F( "fHistEta" , "Pseudo-rap. of casc. candidates ; #eta ; Counts", 120, -3.0, 3.0);
 498         fListHistCascade->Add(fHistEta);
 499 }
 500
 501 if(! fHistTheta ){
 502         fHistTheta  = new TH1F( "fHistTheta" , "#theta of casc. candidates ; #theta (deg) ; Counts", 180, 0., 180.0);
 503         fListHistCascade->Add(fHistTheta);
 504 }
 505
 506 if(! fHistPhi ){
 507         fHistPhi  = new TH1F( "fHistPhi" , "#phi of casc. candidates ; #phi (deg) ; Counts", 360, 0., 360.);
 508         fListHistCascade->Add(fHistPhi);
 509 }
 510
 511
 512 if(! f2dHistArmenteros) {
 513         f2dHistArmenteros = new TH2F( "f2dHistArmenteros", "#alpha_{Arm}(casc. cand.) Vs Pt_{Arm}(casc. cand.); #alpha_{Arm} ; Pt_{Arm} (GeV/c)", 140, -1.2, 1.2, 300, 0., 0.3);
 514         fListHistCascade->Add(f2dHistArmenteros);
 515 }
 516
 517 //-------
 518
 519 if(! f2dHistEffMassLambdaVsEffMassXiMinus) {
 520         f2dHistEffMassLambdaVsEffMassXiMinus = new TH2F( "f2dHistEffMassLambdaVsEffMassXiMinus", "M_{#Lambda} Vs M_{#Xi^{-} candidates} ; Inv. M_{#Lambda^{0}} (GeV/c^{2}) ; M( #Lambda , #pi^{-} ) (GeV/c^{2})", 300, 1.1,1.13, 200, 1.2, 2.0);
 521         fListHistCascade->Add(f2dHistEffMassLambdaVsEffMassXiMinus);
 522 }
 523
 524 if(! f2dHistEffMassXiVsEffMassOmegaMinus) {
 525         f2dHistEffMassXiVsEffMassOmegaMinus = new TH2F( "f2dHistEffMassXiVsEffMassOmegaMinus", "M_{#Xi^{-} candidates} Vs M_{#Omega^{-} candidates} ; M( #Lambda , #pi^{-} ) (GeV/c^{2}) ; M( #Lambda , K^{-} ) (GeV/c^{2})", 200, 1.2, 2.0, 250, 1.5, 2.5);
 526         fListHistCascade->Add(f2dHistEffMassXiVsEffMassOmegaMinus);
 527 }
 528
 529 if(! f2dHistEffMassLambdaVsEffMassXiPlus) {
 530         f2dHistEffMassLambdaVsEffMassXiPlus = new TH2F( "f2dHistEffMassLambdaVsEffMassXiPlus", "M_{#Lambda} Vs M_{#Xi^{+} candidates} ; Inv. M_{#Lambda^{0}} (GeV/c^{2}) ; M( #Lambda , #pi^{+} ) (GeV/c^{2})", 300, 1.1,1.13, 200, 1.2, 2.0);
 531         fListHistCascade->Add(f2dHistEffMassLambdaVsEffMassXiPlus);
 532 }
 533
 534 if(! f2dHistEffMassXiVsEffMassOmegaPlus) {
 535         f2dHistEffMassXiVsEffMassOmegaPlus = new TH2F( "f2dHistEffMassXiVsEffMassOmegaPlus", "M_{#Xi^{+} candidates} Vs M_{#Omega^{+} candidates} ; M( #Lambda , #pi^{+} ) (GeV/c^{2}) ; M( #Lambda , K^{+} ) (GeV/c^{2})", 200, 1.2, 2.0, 250, 1.5, 2.5);
 536         fListHistCascade->Add(f2dHistEffMassXiVsEffMassOmegaPlus);
 537 }
 538
 539 //-------
 540
 541 if(! f2dHistXiRadiusVsEffMassXiMinus) {
 542         f2dHistXiRadiusVsEffMassXiMinus = new TH2F( "f2dHistXiRadiusVsEffMassXiMinus", "Transv. R_{Xi Decay} Vs M_{#Xi^{-} candidates}; r_{cascade} (cm); M( #Lambda , #pi^{-} ) (GeV/c^{2}) ", 450, 0., 45.0, 200, 1.2, 2.0);
 543         fListHistCascade->Add(f2dHistXiRadiusVsEffMassXiMinus);
 544 }
 545
 546 if(! f2dHistXiRadiusVsEffMassXiPlus) {
 547         f2dHistXiRadiusVsEffMassXiPlus = new TH2F( "f2dHistXiRadiusVsEffMassXiPlus", "Transv. R_{Xi Decay} Vs M_{#Xi^{+} candidates}; r_{cascade} (cm); M( #Lambda , #pi^{+} ) (GeV/c^{2}) ", 450, 0., 45.0, 200, 1.2, 2.0);
 548         fListHistCascade->Add(f2dHistXiRadiusVsEffMassXiPlus);
 549 }
 550
 551 if(! f2dHistXiRadiusVsEffMassOmegaMinus) {
 552         f2dHistXiRadiusVsEffMassOmegaMinus = new TH2F( "f2dHistXiRadiusVsEffMassOmegaMinus", "Transv. R_{Xi Decay} Vs M_{#Omega^{-} candidates}; r_{cascade} (cm); M( #Lambda , K^{-} ) (GeV/c^{2}) ", 450, 0., 45.0, 250, 1.5, 2.5);
 553         fListHistCascade->Add(f2dHistXiRadiusVsEffMassOmegaMinus);
 554 }
 555
 556 if(! f2dHistXiRadiusVsEffMassOmegaPlus) {
 557         f2dHistXiRadiusVsEffMassOmegaPlus = new TH2F( "f2dHistXiRadiusVsEffMassOmegaPlus", "Transv. R_{Xi Decay} Vs M_{#Omega^{+} candidates}; r_{cascade} (cm); M( #Lambda , K^{+} ) (GeV/c^{2}) ", 450, 0., 45.0, 250, 1.5, 2.5);
 558         fListHistCascade->Add(f2dHistXiRadiusVsEffMassOmegaPlus);
 559 }
 560
 561
 562 // Part 2 : Raw material for yield extraction -------
 563
 564 if(! f3dHistXiPtVsEffMassVsYXiMinus) {
 565         f3dHistXiPtVsEffMassVsYXiMinus = new TH3F( "f3dHistXiPtVsEffMassVsYXiMinus", "Pt_{cascade} Vs M_{#Xi^{-} candidates} Vs Y_{#Xi}; Pt_{cascade} (GeV/c); M( #Lambda , #pi^{-} ) (GeV/c^{2}) ;Y_{#Xi} ", 100, 0., 10.0, 400, 1.2, 2.0, 48, -1.2,1.2);
 566         fListHistCascade->Add(f3dHistXiPtVsEffMassVsYXiMinus);
 567 }
 568
 569 if(! f3dHistXiPtVsEffMassVsYXiPlus) {
 570         f3dHistXiPtVsEffMassVsYXiPlus = new TH3F( "f3dHistXiPtVsEffMassVsYXiPlus", "Pt_{cascade} Vs M_{#Xi^{+} candidates} Vs Y_{#Xi}; Pt_{cascade} (GeV/c); M( #Lambda , #pi^{+} ) (GeV/c^{2}); Y_{#Xi}", 100, 0., 10.0, 400, 1.2, 2.0, 48, -1.2,1.2);
 571         fListHistCascade->Add(f3dHistXiPtVsEffMassVsYXiPlus);
 572 }
 573
 574 if(! f3dHistXiPtVsEffMassVsYOmegaMinus) {
 575         f3dHistXiPtVsEffMassVsYOmegaMinus = new TH3F( "f3dHistXiPtVsEffMassVsYOmegaMinus", "Pt_{cascade} Vs M_{#Omega^{-} candidates} Vs Y_{#Omega}; Pt_{cascade} (GeV/c); M( #Lambda , K^{-} ) (GeV/c^{2}); Y_{#Omega}", 100, 0., 10.0, 500, 1.5, 2.5, 48, -1.2,1.2);
 576         fListHistCascade->Add(f3dHistXiPtVsEffMassVsYOmegaMinus);
 577 }
 578
 579 if(! f3dHistXiPtVsEffMassVsYOmegaPlus) {
 580         f3dHistXiPtVsEffMassVsYOmegaPlus = new TH3F( "f3dHistXiPtVsEffMassVsYOmegaPlus", "Pt_{cascade} Vs M_{#Omega^{+} candidates} Vs Y_{#Omega}; Pt_{cascade} (GeV/c); M( #Lambda , K^{+} ) (GeV/c^{2}); Y_{#Omega}", 100, 0., 10.0, 500, 1.5, 2.5, 48, -1.2,1.2);
 581         fListHistCascade->Add(f3dHistXiPtVsEffMassVsYOmegaPlus);
 582 }
 583
 584 //--
 585 if(! f3dHistXiPtVsEffMassVsYWithCombPIDXiMinus) {
 586         f3dHistXiPtVsEffMassVsYWithCombPIDXiMinus = new TH3F( "f3dHistXiPtVsEffMassVsYWithCombPIDXiMinus", "Pt_{cascade} Vs M_{#Xi^{-} candidates, with PID} Vs Y_{#Xi}; Pt_{cascade} (GeV/c); M( #Lambda , #pi^{-} ) (GeV/c^{2}) ;Y_{#Xi} ", 100, 0., 10.0, 400, 1.2, 2.0, 48, -1.2,1.2);
 587         fListHistCascade->Add(f3dHistXiPtVsEffMassVsYWithCombPIDXiMinus);
 588 }
 589
 590 if(! f3dHistXiPtVsEffMassVsYWithCombPIDXiPlus) {
 591         f3dHistXiPtVsEffMassVsYWithCombPIDXiPlus = new TH3F( "f3dHistXiPtVsEffMassVsYWithCombPIDXiPlus", "Pt_{cascade} Vs M_{#Xi^{+} candidates, with PID} Vs Y_{#Xi}; Pt_{cascade} (GeV/c); M( #Lambda , #pi^{+} ) (GeV/c^{2}); Y_{#Xi}", 100, 0., 10.0, 400, 1.2, 2.0, 48, -1.2,1.2);
 592         fListHistCascade->Add(f3dHistXiPtVsEffMassVsYWithCombPIDXiPlus);
 593 }
 594
 595 if(! f3dHistXiPtVsEffMassVsYWithCombPIDOmegaMinus) {
 596         f3dHistXiPtVsEffMassVsYWithCombPIDOmegaMinus = new TH3F( "f3dHistXiPtVsEffMassVsYWithCombPIDOmegaMinus", "Pt_{cascade} Vs M_{#Omega^{-} candidates, with PID} Vs Y_{#Omega}; Pt_{cascade} (GeV/c); M( #Lambda , K^{-} ) (GeV/c^{2}); Y_{#Omega}", 100, 0., 10.0, 500, 1.5, 2.5, 48, -1.2,1.2);
 597         fListHistCascade->Add(f3dHistXiPtVsEffMassVsYWithCombPIDOmegaMinus);
 598 }
 599
 600 if(! f3dHistXiPtVsEffMassVsYWithCombPIDOmegaPlus) {
 601         f3dHistXiPtVsEffMassVsYWithCombPIDOmegaPlus = new TH3F( "f3dHistXiPtVsEffMassVsYWithCombPIDOmegaPlus", "Pt_{cascade} Vs M_{#Omega^{+} candidates, with PID} Vs Y_{#Omega}; Pt_{cascade} (GeV/c); M( #Lambda , K^{+} ) (GeV/c^{2}); Y_{#Omega}", 100, 0., 10.0, 500, 1.5, 2.5, 48, -1.2,1.2);
 602         fListHistCascade->Add(f3dHistXiPtVsEffMassVsYWithCombPIDOmegaPlus);
 603 }
 604
 605 //--
 606 if(! f3dHistXiPtVsEffMassVsYWith2CombPIDOmegaMinus) {
 607         f3dHistXiPtVsEffMassVsYWith2CombPIDOmegaMinus = new TH3F( "f3dHistXiPtVsEffMassVsYWith2CombPIDOmegaMinus", "Pt_{cascade} Vs M_{#Omega^{-} candidates, with 2 PID} Vs Y_{#Omega}; Pt_{cascade} (GeV/c); M( #Lambda , K^{-} ) (GeV/c^{2}); Y_{#Omega}", 100, 0., 10.0, 500, 1.5, 2.5, 48, -1.2,1.2);
 608         fListHistCascade->Add(f3dHistXiPtVsEffMassVsYWith2CombPIDOmegaMinus);
 609 }
 610
 611 if(! f3dHistXiPtVsEffMassVsYWith2CombPIDOmegaPlus) {
 612         f3dHistXiPtVsEffMassVsYWith2CombPIDOmegaPlus = new TH3F( "f3dHistXiPtVsEffMassVsYWith2CombPIDOmegaPlus", "Pt_{cascade} Vs M_{#Omega^{+} candidates, with 2 PID} Vs Y_{#Omega}; Pt_{cascade} (GeV/c); M( #Lambda , K^{+} ) (GeV/c^{2}); Y_{#Omega}", 100, 0., 10.0, 500, 1.5, 2.5, 48, -1.2,1.2);
 613         fListHistCascade->Add(f3dHistXiPtVsEffMassVsYWith2CombPIDOmegaPlus);
 614 }
 615
 616
 617 if(! fTpcPidManager){
 618
 619   Double_t lAlephParameters[5] = {0.};
 620         // Reasonable parameters extracted for p-p simulation (LHC09a4) - A.Kalweit
 621         lAlephParameters[0] = 4.23232575531564326e+00;//50*0.76176e-1;
 622         lAlephParameters[1] = 8.68482806165147636e+00;//10.632;
 623         lAlephParameters[2] = 1.34000000000000005e-05;//0.13279e-4;
 624         lAlephParameters[3] = 2.30445734159456084e+00;//1.8631;
 625         lAlephParameters[4] = 2.25624744086878559e+00;//1.9479;
 626
 627   fTpcPidManager = new AliTPCpidESD();
 628   fTpcPidManager->SetBetheBlochParameters(lAlephParameters[0]/50.,
 629                                           lAlephParameters[1],
 630                                           lAlephParameters[2],
 631                                           lAlephParameters[3],
 632                                           lAlephParameters[4]);
 633 }
 634
 635
 636 if(! f3dHistXiPtVsEffMassVsYWithTpcPIDOmegaMinus) {
 637         f3dHistXiPtVsEffMassVsYWithTpcPIDOmegaMinus = new TH3F( "f3dHistXiPtVsEffMassVsYWithTpcPIDOmegaMinus", "Pt_{cascade} Vs M_{#Omega^{-} candidates, with TPC PID} Vs Y_{#Omega}; Pt_{cascade} (GeV/c); M( #Lambda , K^{-} ) (GeV/c^{2}); Y_{#Omega}", 100, 0., 10.0, 500, 1.5, 2.5, 48, -1.2,1.2);
 638         fListHistCascade->Add(f3dHistXiPtVsEffMassVsYWithTpcPIDOmegaMinus);
 639 }
 640
 641
 642
 643 if(! fCFContCascadePIDXiMinus)  {
 644   const Int_t  lNbSteps      =  7 ;
 645   const Int_t  lNbVariables  =  4 ;
 646
 647   //array for the number of bins in each dimension :
 648   Int_t lNbBinsPerVar[4];
 649   lNbBinsPerVar[0] = 100;
 650   lNbBinsPerVar[1] = 400;
 651   lNbBinsPerVar[2] = 48;
 652   lNbBinsPerVar[3] = 200;
 653
 654
 655   fCFContCascadePIDXiMinus = new AliCFContainer("fCFContCascadePIDXiMinus","Pt_{cascade} Vs M_{#Xi^{-} candidates} Vs Y_{#Xi}", lNbSteps, lNbVariables, lNbBinsPerVar );
 656
 657   //setting the bin limits (valid  for v4-18-10-AN)
 658   fCFContCascadePIDXiMinus->SetBinLimits(0,   0.0  ,  10.0 );   // Pt(Cascade)
 659   fCFContCascadePIDXiMinus->SetBinLimits(1,   1.2  ,   2.0 );   // Xi Effective mass
 660   fCFContCascadePIDXiMinus->SetBinLimits(2,  -1.2  ,   1.2 );   // Rapidity
 661   if(fCollidingSystems)
 662         fCFContCascadePIDXiMinus->SetBinLimits(3, 0.0, 20000.0  );    // TrackMultiplicity
 663   else
 664         fCFContCascadePIDXiMinus->SetBinLimits(3, 0.0, 200.0  );     // TrackMultiplicity
 665
 666   // Setting the step title : one per PID case
 667   fCFContCascadePIDXiMinus->SetStepTitle(0, "No PID");
 668   fCFContCascadePIDXiMinus->SetStepTitle(1, "TPC PID / 3-#sigma cut on Bachelor track");
 669   fCFContCascadePIDXiMinus->SetStepTitle(2, "TPC PID / 3-#sigma cut on Bachelor+Baryon tracks");
 670   fCFContCascadePIDXiMinus->SetStepTitle(3, "TPC PID / 3-#sigma cut on Bachelor+Baryon+Meson tracks");
 671   fCFContCascadePIDXiMinus->SetStepTitle(4, "Comb. PID / Bachelor");
 672   fCFContCascadePIDXiMinus->SetStepTitle(5, "Comb. PID / Bachelor+Baryon");
 673   fCFContCascadePIDXiMinus->SetStepTitle(6, "Comb. PID / Bachelor+Baryon+Meson");
 674
 675   // Setting the variable title, per axis
 676   fCFContCascadePIDXiMinus->SetVarTitle(0, "Pt_{cascade} (GeV/c)");
 677   fCFContCascadePIDXiMinus->SetVarTitle(1, "M( #Lambda , #pi^{-} ) (GeV/c^{2})");
 678   fCFContCascadePIDXiMinus->SetVarTitle(2, "Y_{#Xi}");
 679   fCFContCascadePIDXiMinus->SetVarTitle(3, "Track Multiplicity");
 680
 681   fListHistCascade->Add(fCFContCascadePIDXiMinus);
 682
 683 }
 684
 685 if(! fCFContCascadePIDXiPlus)  {
 686   const Int_t  lNbSteps      =  7 ;
 687   const Int_t  lNbVariables  =  4 ;
 688
 689   //array for the number of bins in each dimension :
 690   Int_t lNbBinsPerVar[4];
 691   lNbBinsPerVar[0] = 100;
 692   lNbBinsPerVar[1] = 400;
 693   lNbBinsPerVar[2] = 48;
 694   lNbBinsPerVar[3] = 200;
 695
 696
 697   fCFContCascadePIDXiPlus = new AliCFContainer("fCFContCascadePIDXiPlus","Pt_{cascade} Vs M_{#Xi^{+} candidates} Vs Y_{#Xi}", lNbSteps, lNbVariables, lNbBinsPerVar );
 698
 699
 700   //setting the bin limits (valid  for v4-18-10-AN)
 701   fCFContCascadePIDXiPlus->SetBinLimits(0,   0.0  ,  10.0 );    // Pt(Cascade)
 702   fCFContCascadePIDXiPlus->SetBinLimits(1,   1.2  ,   2.0 );    // Xi Effective mass
 703   fCFContCascadePIDXiPlus->SetBinLimits(2,  -1.2  ,   1.2 );    // Rapidity
 704   if(fCollidingSystems)
 705         fCFContCascadePIDXiPlus->SetBinLimits(3, 0.0, 20000.0  );    // TrackMultiplicity
 706   else
 707         fCFContCascadePIDXiPlus->SetBinLimits(3, 0.0, 200.0  );     // TrackMultiplicity
 708
 709   // Setting the step title : one per PID case
 710   fCFContCascadePIDXiPlus->SetStepTitle(0, "No PID");
 711   fCFContCascadePIDXiPlus->SetStepTitle(1, "TPC PID / 3-#sigma cut on Bachelor track");
 712   fCFContCascadePIDXiPlus->SetStepTitle(2, "TPC PID / 3-#sigma cut on Bachelor+Baryon tracks");
 713   fCFContCascadePIDXiPlus->SetStepTitle(3, "TPC PID / 3-#sigma cut on Bachelor+Baryon+Meson tracks");
 714   fCFContCascadePIDXiPlus->SetStepTitle(4, "Comb. PID / Bachelor");
 715   fCFContCascadePIDXiPlus->SetStepTitle(5, "Comb. PID / Bachelor+Baryon");
 716   fCFContCascadePIDXiPlus->SetStepTitle(6, "Comb. PID / Bachelor+Baryon+Meson");
 717
 718   // Setting the variable title, per axis
 719   fCFContCascadePIDXiPlus->SetVarTitle(0, "Pt_{cascade} (GeV/c)");
 720   fCFContCascadePIDXiPlus->SetVarTitle(1, "M( #Lambda , #pi^{+} ) (GeV/c^{2})");
 721   fCFContCascadePIDXiPlus->SetVarTitle(2, "Y_{#Xi}");
 722   fCFContCascadePIDXiPlus->SetVarTitle(3, "Track Multiplicity");
 723
 724   fListHistCascade->Add(fCFContCascadePIDXiPlus);
 725
 726 }
 727
 728
 729 if(! fCFContCascadePIDOmegaMinus)  {
 730   const Int_t  lNbSteps      =  7 ;
 731   const Int_t  lNbVariables  =  4 ;
 732
 733   //array for the number of bins in each dimension :
 734   Int_t lNbBinsPerVar[4];
 735   lNbBinsPerVar[0] = 100;
 736   lNbBinsPerVar[1] = 500;
 737   lNbBinsPerVar[2] = 48;
 738   lNbBinsPerVar[3] = 200;
 739
 740
 741   fCFContCascadePIDOmegaMinus = new AliCFContainer("fCFContCascadePIDOmegaMinus","Pt_{cascade} Vs M_{#Omega^{-} candidates} Vs Y_{#Omega}", lNbSteps, lNbVariables, lNbBinsPerVar );
 742
 743
 744   //setting the bin limits (valid  for v4-18-10-AN)
 745   fCFContCascadePIDOmegaMinus->SetBinLimits(0,   0.0  ,  10.0 );        // Pt(Cascade)
 746   fCFContCascadePIDOmegaMinus->SetBinLimits(1,   1.5  ,   2.5 );        // Omega Effective mass
 747   fCFContCascadePIDOmegaMinus->SetBinLimits(2,  -1.2  ,   1.2 );        // Rapidity
 748   if(fCollidingSystems)
 749         fCFContCascadePIDOmegaMinus->SetBinLimits(3, 0.0, 20000.0  );    // TrackMultiplicity
 750   else
 751         fCFContCascadePIDOmegaMinus->SetBinLimits(3, 0.0, 200.0  );     // TrackMultiplicity
 752
 753   // Setting the step title : one per PID case
 754   fCFContCascadePIDOmegaMinus->SetStepTitle(0, "No PID");
 755   fCFContCascadePIDOmegaMinus->SetStepTitle(1, "TPC PID / 3-#sigma cut on Bachelor track");
 756   fCFContCascadePIDOmegaMinus->SetStepTitle(2, "TPC PID / 3-#sigma cut on Bachelor+Baryon tracks");
 757   fCFContCascadePIDOmegaMinus->SetStepTitle(3, "TPC PID / 3-#sigma cut on Bachelor+Baryon+Meson tracks");
 758   fCFContCascadePIDOmegaMinus->SetStepTitle(4, "Comb. PID / Bachelor");
 759   fCFContCascadePIDOmegaMinus->SetStepTitle(5, "Comb. PID / Bachelor+Baryon");
 760   fCFContCascadePIDOmegaMinus->SetStepTitle(6, "Comb. PID / Bachelor+Baryon+Meson");
 761
 762   // Setting the variable title, per axis
 763   fCFContCascadePIDOmegaMinus->SetVarTitle(0, "Pt_{cascade} (GeV/c)");
 764   fCFContCascadePIDOmegaMinus->SetVarTitle(1, "M( #Lambda , K^{-} ) (GeV/c^{2})");
 765   fCFContCascadePIDOmegaMinus->SetVarTitle(2, "Y_{#Omega}");
 766   fCFContCascadePIDOmegaMinus->SetVarTitle(3, "Track Multiplicity");
 767
 768   fListHistCascade->Add(fCFContCascadePIDOmegaMinus);
 769
 770 }
 771
 772 if(! fCFContCascadePIDOmegaPlus)  {
 773   const Int_t  lNbSteps      =  7 ;
 774   const Int_t  lNbVariables  =  4 ;
 775
 776   //array for the number of bins in each dimension :
 777   Int_t lNbBinsPerVar[4];
 778   lNbBinsPerVar[0] = 100;
 779   lNbBinsPerVar[1] = 500;
 780   lNbBinsPerVar[2] = 48;
 781   lNbBinsPerVar[3] = 200;
 782
 783
 784   fCFContCascadePIDOmegaPlus = new AliCFContainer("fCFContCascadePIDOmegaPlus","Pt_{cascade} Vs M_{#Omega^{+} candidates} Vs Y_{#Omega}", lNbSteps, lNbVariables, lNbBinsPerVar );
 785
 786
 787   //setting the bin limits (valid  for v4-18-10-AN)
 788   fCFContCascadePIDOmegaPlus->SetBinLimits(0,   0.0  ,  10.0 ); // Pt(Cascade)
 789   fCFContCascadePIDOmegaPlus->SetBinLimits(1,   1.5  ,   2.5 ); // Omega Effective mass
 790   fCFContCascadePIDOmegaPlus->SetBinLimits(2,  -1.2  ,   1.2 ); // Rapidity
 791   if(fCollidingSystems)
 792         fCFContCascadePIDOmegaPlus->SetBinLimits(3, 0.0, 20000.0  );    // TrackMultiplicity
 793   else
 794         fCFContCascadePIDOmegaPlus->SetBinLimits(3, 0.0, 200.0  );     // TrackMultiplicity
 795
 796   // Setting the step title : one per PID case
 797   fCFContCascadePIDOmegaPlus->SetStepTitle(0, "No PID");
 798   fCFContCascadePIDOmegaPlus->SetStepTitle(1, "TPC PID / 3-#sigma cut on Bachelor track");
 799   fCFContCascadePIDOmegaPlus->SetStepTitle(2, "TPC PID / 3-#sigma cut on Bachelor+Baryon tracks");
 800   fCFContCascadePIDOmegaPlus->SetStepTitle(3, "TPC PID / 3-#sigma cut on Bachelor+Baryon+Meson tracks");
 801   fCFContCascadePIDOmegaPlus->SetStepTitle(4, "Comb. PID / Bachelor");
 802   fCFContCascadePIDOmegaPlus->SetStepTitle(5, "Comb. PID / Bachelor+Baryon");
 803   fCFContCascadePIDOmegaPlus->SetStepTitle(6, "Comb. PID / Bachelor+Baryon+Meson");
 804
 805   // Setting the variable title, per axis
 806   fCFContCascadePIDOmegaPlus->SetVarTitle(0, "Pt_{cascade} (GeV/c)");
 807   fCFContCascadePIDOmegaPlus->SetVarTitle(1, "M( #Lambda , K^{+} ) (GeV/c^{2})");
 808   fCFContCascadePIDOmegaPlus->SetVarTitle(2, "Y_{#Omega}");
 809   fCFContCascadePIDOmegaPlus->SetVarTitle(3, "Track Multiplicity");
 810
 811   fListHistCascade->Add(fCFContCascadePIDOmegaPlus);
 812
 813 }
 814
 815
 816
 817
 818
 819
 820
 821
 822
 823
 824
 825 // Part 3 : Towards the optimisation of topological selections -------
 826 if(! fCFContCascadeCuts){
 827
 828         // Container meant to store all the relevant distributions corresponding to the cut variables.
 829         // So far, 19 variables have been identified.
 830         // The following will be done in quite an inelegant way.
 831         // Improvement expected later.
 832   const Int_t  lNbSteps      =  2 ;
 833   const Int_t  lNbVariables  =  18 ;
 834
 835   //array for the number of bins in each dimension :
 836   Int_t lNbBinsPerVar[18];
 837   lNbBinsPerVar[0]  = 25;
 838   lNbBinsPerVar[1]  = 25;
 839   lNbBinsPerVar[2]  = 20;
 840   lNbBinsPerVar[3]  = 40;
 841   lNbBinsPerVar[4]  = 50;
 842   lNbBinsPerVar[5]  = 12;
 843
 844   lNbBinsPerVar[6]  = 20;
 845   lNbBinsPerVar[7]  = 40;
 846   lNbBinsPerVar[8]  = 40;
 847   lNbBinsPerVar[9]  = 25;
 848   lNbBinsPerVar[10] = 25;
 849
 850   lNbBinsPerVar[11] = 60;
 851   lNbBinsPerVar[12] = 50;
 852
 853   lNbBinsPerVar[13] = 20;
 854   lNbBinsPerVar[14] = 20;
 855
 856   lNbBinsPerVar[15] = 40;
 857   lNbBinsPerVar[16] = 40;
 858   lNbBinsPerVar[17] = 35;
 859
 860  fCFContCascadeCuts = new AliCFContainer("fCFContCascadeCuts","Container for Cascade cuts", lNbSteps, lNbVariables, lNbBinsPerVar );
 861
 862
 863   //setting the bin limits (valid for v4-18-10-AN on)
 864   fCFContCascadeCuts->SetBinLimits(0,    0.0  ,  0.25 );        // DcaXiDaughters
 865   fCFContCascadeCuts->SetBinLimits(1,    0.0  ,  0.25 );        // DcaBachToPrimVertexXi
 866   fCFContCascadeCuts->SetBinLimits(2,    0.995,  1.0  );        // XiCosineOfPointingAngle
 867   fCFContCascadeCuts->SetBinLimits(3,    0.0  ,  4.0  );        // XiRadius
 868   fCFContCascadeCuts->SetBinLimits(4,    1.1  ,  1.15  );       // InvMassLambdaAsCascDghter
 869   fCFContCascadeCuts->SetBinLimits(5,    0.0  ,  0.6  );        // DcaV0DaughtersXi
 870   fCFContCascadeCuts->SetBinLimits(6,    0.98 ,  1.0  );        // V0CosineOfPointingAngleXi
 871   fCFContCascadeCuts->SetBinLimits(7,    0.0  , 20.0  );        // V0RadiusXi
 872   fCFContCascadeCuts->SetBinLimits(8,    0.0  ,  1.0  );        // DcaV0ToPrimVertexXi
 873   fCFContCascadeCuts->SetBinLimits(9,    0.0  ,  2.5  );        // DcaPosToPrimVertexXi
 874   fCFContCascadeCuts->SetBinLimits(10,   0.0  ,  2.5  );        // DcaNegToPrimVertexXi
 875   fCFContCascadeCuts->SetBinLimits(11,   1.25 ,  1.45  );       // InvMassXi
 876   fCFContCascadeCuts->SetBinLimits(12,   1.6  ,  1.8  );        // InvMassOmega
 877   fCFContCascadeCuts->SetBinLimits(13,   0.0  , 10.0  );        // XiTransvMom
 878   fCFContCascadeCuts->SetBinLimits(14, -10.0  , 10.0  );        // BestPrimaryVtxPosZ
 879   if(fCollidingSystems){
 880         fCFContCascadeCuts->SetBinLimits(15,   0.0, 10000.0  );    // TrackMultiplicity
 881         fCFContCascadeCuts->SetBinLimits(16,   0.0, 10000.0  );    // SPDTrackletsMultiplicity
 882   }
 883   else{
 884         fCFContCascadeCuts->SetBinLimits(15,   0.0, 200.0  );     // TrackMultiplicity
 885         fCFContCascadeCuts->SetBinLimits(16,   0.0, 200.0  );     // SPDTrackletsMultiplicity
 886   }
 887   fCFContCascadeCuts->SetBinLimits(17,  25.0  ,165.0  );        // BachTPCClusters
 888
 889
 890
 891   // Setting the number of steps : one for negative cascades (Xi- and Omega-), another for the positve cascades(Xi+ and Omega+)
 892   fCFContCascadeCuts->SetStepTitle(0, "Negative Cascades");
 893   fCFContCascadeCuts->SetStepTitle(1, "Positive Cascades");
 894
 895   // Setting the variable title, per axis
 896   // fCFContCascadeCuts->SetVarTitle(0,  "Chi2Xi");
 897   fCFContCascadeCuts->SetVarTitle(0,  "DcaXiDaughters");
 898   fCFContCascadeCuts->SetVarTitle(1,  "DcaBachToPrimVertexXi");
 899   fCFContCascadeCuts->SetVarTitle(2,  "XiCosineOfPointingAngle");
 900   fCFContCascadeCuts->SetVarTitle(3,  "XiRadius");
 901   fCFContCascadeCuts->SetVarTitle(4,  "InvMassLambdaAsCascDghter");
 902    // fCFContCascadeCuts->SetVarTitle(0,  "V0Chi2Xi");
 903   fCFContCascadeCuts->SetVarTitle(5,  "DcaV0DaughtersXi");
 904
 905   fCFContCascadeCuts->SetVarTitle(6,  "V0CosineOfPointingAngleXi");
 906   fCFContCascadeCuts->SetVarTitle(7,  "V0RadiusXi");
 907   fCFContCascadeCuts->SetVarTitle(8,  "DcaV0ToPrimVertexXi");
 908   fCFContCascadeCuts->SetVarTitle(9,  "DcaPosToPrimVertexXi");
 909   fCFContCascadeCuts->SetVarTitle(10, "DcaNegToPrimVertexXi");
 910
 911   fCFContCascadeCuts->SetVarTitle(11, "InvMassXi");
 912   fCFContCascadeCuts->SetVarTitle(12, "InvMassOmega");
 913
 914   fCFContCascadeCuts->SetVarTitle(13, "XiTransvMom");
 915   //fCFContCascadeCuts->SetVarTitle(14, "V0toXiCosineOfPointingAngle");
 916   fCFContCascadeCuts->SetVarTitle(14, "BestPrimaryVtxPosZ");
 917
 918   fCFContCascadeCuts->SetVarTitle(15, "TrackMultiplicity");
 919   fCFContCascadeCuts->SetVarTitle(16, "SPDTrackletsMultiplicity");
 920   fCFContCascadeCuts->SetVarTitle(17, "BachTPCClusters");
 921
 922   fListHistCascade->Add(fCFContCascadeCuts);
 923 }
 924
 925
 926
 927 // Part 4 : Angular correlation study -------
 928
 929 if(! fHnSpAngularCorrXiMinus){
 930         // Delta Phi(Casc,any trck) Vs Delta Eta(Casc,any trck) Vs Casc Pt Vs Pt of the tracks
 931         // Delta Phi  = 360 bins de -180., 180.
 932         // Delta Eta  = 120 bins de -3.0, 3.0
 933         // Pt Cascade = 100 bins de 0., 10.0,
 934         // Pt track = 150 bins de 0., 15.0
 935
 936    Int_t    bins[5] = { 360, 120, 100, 150, 40};
 937    Double_t xmin[5] = {-50., -3.,  0.,  0., 1.30};
 938    Double_t xmax[5] = { 310., 3., 10., 15., 1.34};
 939    fHnSpAngularCorrXiMinus = new THnSparseF("fHnSpAngularCorrXiMinus", "Angular Correlation for #Xi^{-}:", 5, bins, xmin, xmax);
 940         fHnSpAngularCorrXiMinus->GetAxis(0)->SetTitle(" #Delta#phi(Casc,Track) (deg)");
 941         fHnSpAngularCorrXiMinus->GetAxis(1)->SetTitle(" #Delta#eta(Casc,Track)");
 942         fHnSpAngularCorrXiMinus->GetAxis(2)->SetTitle(" Pt_{Casc} (GeV/c)");
 943         fHnSpAngularCorrXiMinus->GetAxis(3)->SetTitle(" Pt_{any track} (GeV/c)");
 944         fHnSpAngularCorrXiMinus->GetAxis(4)->SetTitle(" Eff. Inv Mass (GeV/c^{2})");
 945         fHnSpAngularCorrXiMinus->Sumw2();
 946    fListHistCascade->Add(fHnSpAngularCorrXiMinus);
 947 }
 948
 949 if(! fHnSpAngularCorrXiPlus){
 950         // Delta Phi(Casc,any trck) Vs Delta Eta(Casc,any trck) Vs Casc Pt Vs Pt of the tracks
 951         // Delta Phi  = 360 bins de -180., 180.
 952         // Delta Eta  = 120 bins de -3.0, 3.0
 953         // Pt Cascade = 100 bins de 0., 10.0,
 954         // Pt track = 150 bins de 0., 15.0
 955    Int_t    bins[5] = { 360, 120, 100, 150, 40};
 956    Double_t xmin[5] = {-50., -3.,  0.,  0., 1.30};
 957    Double_t xmax[5] = { 310., 3., 10., 15., 1.34};
 958    fHnSpAngularCorrXiPlus = new THnSparseF("fHnSpAngularCorrXiPlus", "Angular Correlation for #Xi^{+}:", 5, bins, xmin, xmax);
 959         fHnSpAngularCorrXiPlus->GetAxis(0)->SetTitle(" #Delta#phi(Casc,Track) (deg)");
 960         fHnSpAngularCorrXiPlus->GetAxis(1)->SetTitle(" #Delta#eta(Casc,Track)");
 961         fHnSpAngularCorrXiPlus->GetAxis(2)->SetTitle(" Pt_{Casc} (GeV/c)");
 962         fHnSpAngularCorrXiPlus->GetAxis(3)->SetTitle(" Pt_{any track} (GeV/c)");
 963         fHnSpAngularCorrXiPlus->GetAxis(4)->SetTitle(" Eff. Inv Mass (GeV/c^{2})");
 964         fHnSpAngularCorrXiPlus->Sumw2();
 965    fListHistCascade->Add(fHnSpAngularCorrXiPlus);
 966 }
 967
 968 if(! fHnSpAngularCorrOmegaMinus){
 969         // Delta Phi(Casc,any trck) Vs Delta Eta(Casc,any trck) Vs Casc Pt Vs Pt of the tracks
 970         // Delta Phi  = 360 bins de -180., 180.
 971         // Delta Eta  = 120 bins de -3.0, 3.0
 972         // Pt Cascade = 100 bins de 0., 10.0,
 973         // Pt track = 150 bins de 0., 15.0
 974
 975    Int_t    bins[5] = { 360, 120, 100, 150, 40};
 976    Double_t xmin[5] = {-50., -3.,  0.,  0., 1.65};
 977    Double_t xmax[5] = { 310., 3., 10., 15., 1.69};
 978    fHnSpAngularCorrOmegaMinus = new THnSparseF("fHnSpAngularCorrOmegaMinus", "Angular Correlation for #Omega^{-}:", 5, bins, xmin, xmax);
 979         fHnSpAngularCorrOmegaMinus->GetAxis(0)->SetTitle(" #Delta#phi(Casc,Track) (deg)");
 980         fHnSpAngularCorrOmegaMinus->GetAxis(1)->SetTitle(" #Delta#eta(Casc,Track)");
 981         fHnSpAngularCorrOmegaMinus->GetAxis(2)->SetTitle(" Pt_{Casc} (GeV/c)");
 982         fHnSpAngularCorrOmegaMinus->GetAxis(3)->SetTitle(" Pt_{any track} (GeV/c)");
 983         fHnSpAngularCorrOmegaMinus->GetAxis(4)->SetTitle(" Eff. Inv Mass (GeV/c^{2})");
 984         fHnSpAngularCorrOmegaMinus->Sumw2();
 985    fListHistCascade->Add(fHnSpAngularCorrOmegaMinus);
 986 }
 987
 988 if(! fHnSpAngularCorrOmegaPlus){
 989         // Delta Phi(Casc,any trck) Vs Delta Eta(Casc,any trck) Vs Casc Pt Vs Pt of the tracks
 990         // Delta Phi  = 360 bins de -180., 180.
 991         // Delta Eta  = 120 bins de -3.0, 3.0
 992         // Pt Cascade = 100 bins de 0., 10.0,
 993         // Pt track = 150 bins de 0., 15.0
 994    Int_t    bins[5] = { 360, 120, 100, 150, 40};
 995    Double_t xmin[5] = {-50., -3.,  0.,  0., 1.65};
 996    Double_t xmax[5] = { 310., 3., 10., 15., 1.69};
 997    fHnSpAngularCorrOmegaPlus = new THnSparseF("fHnSpAngularCorrOmegaPlus", "Angular Correlation for #Omega^{+}:", 5, bins, xmin, xmax);
 998         fHnSpAngularCorrOmegaPlus->GetAxis(0)->SetTitle(" #Delta#phi(Casc,Track) (deg)");
 999         fHnSpAngularCorrOmegaPlus->GetAxis(1)->SetTitle(" #Delta#eta(Casc,Track)");
1000         fHnSpAngularCorrOmegaPlus->GetAxis(2)->SetTitle(" Pt_{Casc} (GeV/c)");
1001         fHnSpAngularCorrOmegaPlus->GetAxis(3)->SetTitle(" Pt_{any track} (GeV/c)");
1002         fHnSpAngularCorrOmegaPlus->GetAxis(4)->SetTitle(" Eff. Inv Mass (GeV/c^{2})");
1003         fHnSpAngularCorrOmegaPlus->Sumw2();
1004    fListHistCascade->Add(fHnSpAngularCorrOmegaPlus);
1005 }
1006
1007
1008 }// end UserCreateOutputObjects
1009
1010
1011
1012
1013
1014
1015 //________________________________________________________________________
1016 void AliAnalysisTaskCheckCascade::UserExec(Option_t *)
1017 {
1018   // Main loop
1019   // Called for each event
1020
1021         AliESDEvent *lESDevent = 0x0;
1022         AliAODEvent *lAODevent = 0x0;
1023         Int_t ncascades          = -1;
1024         Int_t nTrackMultiplicity = -1;
1025
1026
1027   // Connect to the InputEvent
1028   // After these lines, we should have an ESD/AOD event + the number of cascades in it.
1029
1030   if(fAnalysisType == "ESD"){
1031         lESDevent = dynamic_cast<AliESDEvent*>( InputEvent() );
1032         if (!lESDevent) {
1033                 Printf("ERROR: lESDevent not available \n");
1034                 return;
1035         }
1036         ncascades = lESDevent->GetNumberOfCascades();
1037   }
1038
1039   if(fAnalysisType == "AOD"){
1040           lAODevent = dynamic_cast<AliAODEvent*>( InputEvent() );
1041         if (!lAODevent) {
1042                 Printf("ERROR: lAODevent not available \n");
1043                 return;
1044         }
1045         ncascades = lAODevent->GetNumberOfCascades();
1046         // Printf("Number of cascade(s) = %d \n", ncascades);
1047   }
1048
1049         nTrackMultiplicity = (InputEvent())->GetNumberOfTracks();
1050
1051
1052   //-------------------------------------------------
1053   // O - Cascade vertexer (ESD)
1054
1055         //   if(fAnalysisType == "ESD" ){
1056         //      lESDevent->ResetCascades();
1057         //      AliCascadeVertexer CascVtxer;
1058         //      CascVtxer.V0sTracks2CascadeVertices(lESDevent);
1059         //     }
1060
1061
1062   // ---------------------------------------------------------------
1063   // I - General histos (filled for any event)
1064
1065         fHistTrackMultiplicity  ->Fill( nTrackMultiplicity );
1066         fHistCascadeMultiplicity->Fill( ncascades );
1067
1068
1069   // ---------------------------------------------------------------
1070   // II - Calcultaion Part dedicated to Xi vertices
1071
1072   for (Int_t iXi = 0; iXi < ncascades; iXi++)
1073   {// This is the begining of the Cascade loop (ESD or AOD)
1074
1075     // -------------------------------------
1076     // II.Init - Initialisation of the local variables that will be needed for ESD/AOD
1077
1078
1079         // - 0th part of initialisation : around primary vertex ...
1080         Short_t  lStatusTrackingPrimVtx  = -2;
1081         Double_t lTrkgPrimaryVtxPos[3]   = {-100.0, -100.0, -100.0};
1082         Double_t lTrkgPrimaryVtxRadius3D = -500.0;
1083
1084         Double_t lBestPrimaryVtxPos[3]   = {-100.0, -100.0, -100.0};
1085         Double_t lBestPrimaryVtxRadius3D = -500.0;
1086
1087         // - 1st part of initialisation : variables needed to store AliESDCascade data members
1088         Double_t lEffMassXi      = 0. ;
1089         Double_t lChi2Xi         = 0. ;
1090         Double_t lDcaXiDaughters = 0. ;
1091         Double_t lXiCosineOfPointingAngle = 0. ;
1092         Double_t lPosXi[3] = { -1000.0, -1000.0, -1000.0 };
1093         Double_t lXiRadius = 0. ;
1094
1095
1096         // - 2nd part of initialisation : about V0 part in cascades
1097         Double_t lInvMassLambdaAsCascDghter = 0.;
1098         Double_t lV0Chi2Xi         = 0. ;
1099         Double_t lDcaV0DaughtersXi = 0.;
1100
1101         Double_t lDcaBachToPrimVertexXi = 0., lDcaV0ToPrimVertexXi = 0.;
1102         Double_t lDcaPosToPrimVertexXi  = 0.;
1103         Double_t lDcaNegToPrimVertexXi  = 0.;
1104         Double_t lV0CosineOfPointingAngleXi = 0. ;
1105         Double_t lPosV0Xi[3] = { -1000. , -1000., -1000. }; // Position of VO coming from cascade
1106         Double_t lV0RadiusXi = -1000.0;
1107         Double_t lV0quality  = 0.;
1108
1109
1110         // - 3rd part of initialisation : Effective masses
1111         Double_t lInvMassXiMinus    = 0.;
1112         Double_t lInvMassXiPlus     = 0.;
1113         Double_t lInvMassOmegaMinus = 0.;
1114         Double_t lInvMassOmegaPlus  = 0.;
1115
1116         // - 4th part of initialisation : PID treatment
1117         Bool_t   lIsPosInXiProton      = kFALSE;
1118         Bool_t   lIsPosInXiPion        = kFALSE;
1119         Bool_t   lIsPosInOmegaProton   = kFALSE;
1120         Bool_t   lIsPosInOmegaPion     = kFALSE;
1121
1122         Bool_t   lIsNegInXiProton      = kFALSE;
1123         Bool_t   lIsNegInXiPion        = kFALSE;
1124         Bool_t   lIsNegInOmegaProton   = kFALSE;
1125         Bool_t   lIsNegInOmegaPion     = kFALSE;
1126
1127         Bool_t   lIsBachelorKaon       = kFALSE;
1128         Bool_t   lIsBachelorPion       = kFALSE;
1129
1130         Bool_t   lIsBachelorKaonForTPC = kFALSE; // For ESD only ...//FIXME : wait for availability in AOD
1131         Bool_t   lIsBachelorPionForTPC = kFALSE; // For ESD only ...
1132         Bool_t   lIsNegPionForTPC      = kFALSE; // For ESD only ...
1133         Bool_t   lIsPosPionForTPC      = kFALSE; // For ESD only ...
1134         Bool_t   lIsNegProtonForTPC    = kFALSE; // For ESD only ...
1135         Bool_t   lIsPosProtonForTPC    = kFALSE; // For ESD only ...
1136
1137         // - 5th part of initialisation : extra info for QA
1138         Double_t lXiMomX       = 0., lXiMomY = 0., lXiMomZ = 0.;
1139         Double_t lXiTransvMom  = 0. ;
1140         Double_t lXiTotMom     = 0. ;
1141
1142         Double_t lBachMomX       = 0., lBachMomY  = 0., lBachMomZ   = 0.;
1143         Double_t lBachTransvMom  = 0.;
1144         Double_t lBachTotMom     = 0.;
1145
1146         Short_t  lChargeXi = -2;
1147         Double_t lV0toXiCosineOfPointingAngle = 0. ;
1148
1149         Double_t lRapXi   = -20.0, lRapOmega = -20.0,  lEta = -20.0, lTheta = 360., lPhi = 720. ;
1150         Double_t lAlphaXi = -200., lPtArmXi  = -200.0;
1151
1152         // - 6th part of initialisation : variables for the AliCFContainer dedicated to cascade cut optmisiation
1153         Int_t    lSPDTrackletsMultiplicity = -1;
1154         Int_t    lBachTPCClusters          = -1;
1155
1156         // - 7th part of initialisation : variables needed for Angular correlations
1157         TVector3 lTVect3MomXi(0.,0.,0.);
1158         Int_t    lArrTrackID[3] = {-1, -1, -1};
1159
1160
1161   if(fAnalysisType == "ESD"){
1162
1163   // -------------------------------------
1164   // II.ESD - Calcultaion Part dedicated to Xi vertices (ESD)
1165
1166         AliESDcascade *xi = lESDevent->GetCascade(iXi);
1167         if (!xi) continue;
1168
1169         // Just to know which file is currently open : locate the file containing Xi
1170         // cout << "Name of the file containing Xi candidate(s) :"
1171         //      << fInputHandler->GetTree()->GetCurrentFile()->GetName()
1172         //      << endl;
1173
1174                 // - II.Step 1 : Characteristics of the event : prim. Vtx + magnetic field (ESD)
1175                 //-------------
1176
1177         // For new code (v4-16-Release-Rev06 or trunk)
1178
1179         const AliESDVertex *lPrimaryTrackingVtx = lESDevent->GetPrimaryVertexTracks();
1180                 // get the vtx stored in ESD found with tracks
1181
1182                 lPrimaryTrackingVtx->GetXYZ( lTrkgPrimaryVtxPos );
1183                 lTrkgPrimaryVtxRadius3D = TMath::Sqrt(  lTrkgPrimaryVtxPos[0] * lTrkgPrimaryVtxPos[0] +
1184                                                         lTrkgPrimaryVtxPos[1] * lTrkgPrimaryVtxPos[1] +
1185                                                         lTrkgPrimaryVtxPos[2] * lTrkgPrimaryVtxPos[2] );
1186
1187                 lStatusTrackingPrimVtx = lPrimaryTrackingVtx->GetStatus();
1188
1189
1190         const AliESDVertex *lPrimaryBestVtx = lESDevent->GetPrimaryVertex();
1191                 // get the best primary vertex available for the event
1192                 // As done in AliCascadeVertexer, we keep the one which is the best one available.
1193                 // between : Tracking vertex > SPD vertex > TPC vertex > default SPD vertex
1194                 // This one will be used for next calculations (DCA essentially)
1195
1196                 lPrimaryBestVtx->GetXYZ( lBestPrimaryVtxPos );
1197                 lBestPrimaryVtxRadius3D = TMath::Sqrt(  lBestPrimaryVtxPos[0] * lBestPrimaryVtxPos[0] +
1198                                                         lBestPrimaryVtxPos[1] * lBestPrimaryVtxPos[1] +
1199                                                         lBestPrimaryVtxPos[2] * lBestPrimaryVtxPos[2] );
1200
1201
1202         // For older evts
1203
1204         //      const AliESDVertex *lPrimaryTrackingVtx = lESDevent->GetPrimaryVertexTracks();
1205         //              get the vtx stored in ESD found with tracks
1206         //      Double_t lTrkgPrimaryVtxPos[3] = {-100.0, -100.0, -100.0};
1207         //              lPrimaryTrackingVtx->GetXYZ( lTrkgPrimaryVtxPos );
1208         //              Double_t lTrkgPrimaryVtxRadius3D = TMath::Sqrt( lTrkgPrimaryVtxPos[0]*lTrkgPrimaryVtxPos[0] +
1209         //                                              lTrkgPrimaryVtxPos[1] * lTrkgPrimaryVtxPos[1] +
1210         //                                              lTrkgPrimaryVtxPos[2] * lTrkgPrimaryVtxPos[2] );
1211         //
1212         //      const AliESDVertex *lPrimarySPDVtx = lESDevent->GetPrimaryVertexSPD();
1213         //              get the vtx found by exclusive use of SPD
1214         //      Double_t lSPDPrimaryVtxPos[3] = {-100.0, -100.0, -100.0};
1215         //              lPrimarySPDVtx->GetXYZ( lSPDPrimaryVtxPos );
1216         //
1217         //      // As done in AliCascadeVertexer, we keep, between both retrieved vertices,
1218         //      // the one which is the best one available.
1219         //      // This one will be used for next calculations (DCA essentially)
1220         //
1221         //              Double_t lBestPrimaryVtxPos[3] = {-100.0, -100.0, -100.0};
1222         //              if( lPrimaryTrackingVtx->GetStatus() ) { // if tracking vtx = ok
1223         //                      lBestPrimaryVtxPos[0] = lTrkgPrimaryVtxPos[0];
1224         //                      lBestPrimaryVtxPos[1] = lTrkgPrimaryVtxPos[1];
1225         //                      lBestPrimaryVtxPos[2] = lTrkgPrimaryVtxPos[2];
1226         //              }
1227         //              else{
1228         //                      lBestPrimaryVtxPos[0] = lSPDPrimaryVtxPos[0];
1229         //                      lBestPrimaryVtxPos[1] = lSPDPrimaryVtxPos[1];
1230         //                      lBestPrimaryVtxPos[2] = lSPDPrimaryVtxPos[2];
1231         //              }
1232         //
1233         //      Double_t lBestPrimaryVtxRadius3D = TMath::Sqrt( lBestPrimaryVtxPos[0]*lBestPrimaryVtxPos[0] +
1234         //                                              lBestPrimaryVtxPos[1] * lBestPrimaryVtxPos[1] +
1235         //                                              lBestPrimaryVtxPos[2] * lBestPrimaryVtxPos[2] );
1236
1237         // Back to normal
1238
1239         Double_t lMagneticField = lESDevent->GetMagneticField( );
1240
1241
1242
1243                 // - II.Step 2 : Assigning the necessary variables for specific AliESDcascade data members (ESD)
1244                 //-------------
1245         lV0quality = 0.;
1246         xi->ChangeMassHypothesis(lV0quality , 3312); // default working hypothesis : cascade = Xi- decay
1247
1248         lEffMassXi                      = xi->GetEffMassXi();
1249         lChi2Xi                         = xi->GetChi2Xi();
1250         lDcaXiDaughters                 = xi->GetDcaXiDaughters();
1251         lXiCosineOfPointingAngle        = xi->GetCascadeCosineOfPointingAngle( lBestPrimaryVtxPos[0],
1252                                                                                 lBestPrimaryVtxPos[1],
1253                                                                                 lBestPrimaryVtxPos[2] );
1254                 // Take care : the best available vertex should be used (like in AliCascadeVertexer)
1255
1256         xi->GetXYZcascade( lPosXi[0],  lPosXi[1], lPosXi[2] );
1257         lXiRadius                       = TMath::Sqrt( lPosXi[0]*lPosXi[0]  +  lPosXi[1]*lPosXi[1] );
1258
1259
1260
1261                 // - II.Step 3 : around the tracks : Bach + V0 (ESD)
1262                 // ~ Necessary variables for ESDcascade data members coming from the ESDv0 part (inheritance)
1263                 //-------------
1264
1265                 UInt_t lIdxPosXi        = (UInt_t) TMath::Abs( xi->GetPindex() );
1266                 UInt_t lIdxNegXi        = (UInt_t) TMath::Abs( xi->GetNindex() );
1267                 UInt_t lBachIdx         = (UInt_t) TMath::Abs( xi->GetBindex() );
1268                         // Care track label can be negative in MC production (linked with the track quality)
1269                         // However = normally, not the case for track index ...
1270
1271         AliESDtrack *pTrackXi           = lESDevent->GetTrack( lIdxPosXi );
1272         AliESDtrack *nTrackXi           = lESDevent->GetTrack( lIdxNegXi );
1273         AliESDtrack *bachTrackXi        = lESDevent->GetTrack( lBachIdx );
1274         if (!pTrackXi || !nTrackXi || !bachTrackXi ) {
1275                 Printf("ERROR: Could not retrieve one of the 3 ESD daughter tracks of the cascade ...");
1276                 continue;
1277         }
1278
1279         lInvMassLambdaAsCascDghter      = xi->GetEffMass();
1280                 // This value shouldn't change, whatever the working hyp. is : Xi-, Xi+, Omega-, Omega+
1281         lDcaV0DaughtersXi               = xi->GetDcaV0Daughters();
1282         lV0Chi2Xi                       = xi->GetChi2V0();
1283
1284         lV0CosineOfPointingAngleXi      = xi->GetV0CosineOfPointingAngle( lBestPrimaryVtxPos[0],
1285                                                                           lBestPrimaryVtxPos[1],
1286                                                                           lBestPrimaryVtxPos[2] );
1287
1288         lDcaV0ToPrimVertexXi            = xi->GetD( lBestPrimaryVtxPos[0],
1289                                                     lBestPrimaryVtxPos[1],
1290                                                     lBestPrimaryVtxPos[2] );
1291
1292         lDcaBachToPrimVertexXi = TMath::Abs( bachTrackXi->GetD( lBestPrimaryVtxPos[0],
1293                                                                 lBestPrimaryVtxPos[1],
1294                                                                 lMagneticField  ) );
1295                                         // Note : AliExternalTrackParam::GetD returns an algebraic value ...
1296
1297                 xi->GetXYZ( lPosV0Xi[0],  lPosV0Xi[1], lPosV0Xi[2] );
1298         lV0RadiusXi             = TMath::Sqrt( lPosV0Xi[0]*lPosV0Xi[0]  +  lPosV0Xi[1]*lPosV0Xi[1] );
1299
1300         lDcaPosToPrimVertexXi   = TMath::Abs( pTrackXi  ->GetD( lBestPrimaryVtxPos[0],
1301                                                                 lBestPrimaryVtxPos[1],
1302                                                                 lMagneticField  )     );
1303
1304         lDcaNegToPrimVertexXi   = TMath::Abs( nTrackXi  ->GetD( lBestPrimaryVtxPos[0],
1305                                                                 lBestPrimaryVtxPos[1],
1306                                                                 lMagneticField  )     );
1307
1308
1309                 // - II.Step 3' : extra-selection for cascade candidates
1310                 // Towards optimisation of AA selection
1311         Bool_t kExtraSelections = kFALSE;
1312
1313         if(kExtraSelections){
1314                 // if(lChi2Xi > 2000) continue;
1315                 // if(lV0Chi2Xi > 2000) continue;
1316
1317                 if(lDcaXiDaughters > 0.05) continue; // > 0.1 by default
1318                 //if(lXiCosineOfPointingAngle < 0.999 ) continue;
1319                 if(lXiRadius < 1.0) continue;
1320                 if(lXiRadius > 100) continue;
1321                 if(TMath::Abs(lInvMassLambdaAsCascDghter-1.11568) > 0.007) continue;
1322                 if(lDcaV0DaughtersXi > 0.3) continue;
1323
1324                 if(lV0CosineOfPointingAngleXi > 0.9999) continue;
1325                 //if(lDcaV0ToPrimVertexXi < 0.09) continue;
1326                 if(lDcaBachToPrimVertexXi < 0.04) continue;
1327
1328                 if(lV0RadiusXi < 1.0) continue;
1329                 if(lV0RadiusXi > 100) continue;
1330                 //if(lDcaPosToPrimVertexXi < 0.6) continue;
1331                 //if(lDcaNegToPrimVertexXi < 0.6) continue;
1332         }
1333
1334
1335
1336                 // - II.Step 4 : around effective masses (ESD)
1337                 // ~ change mass hypotheses to cover all the possibilities :  Xi-/+, Omega -/+
1338                 //-------------
1339
1340
1341         if( bachTrackXi->Charge() < 0 ) {
1342                 lV0quality = 0.;
1343                 xi->ChangeMassHypothesis(lV0quality , 3312);
1344                         // Calculate the effective mass of the Xi- candidate.
1345                         // pdg code 3312 = Xi-
1346                 lInvMassXiMinus = xi->GetEffMassXi();
1347
1348                 lV0quality = 0.;
1349                 xi->ChangeMassHypothesis(lV0quality , 3334);
1350                         // Calculate the effective mass of the Xi- candidate.
1351                         // pdg code 3334 = Omega-
1352                 lInvMassOmegaMinus = xi->GetEffMassXi();
1353
1354                 lV0quality = 0.;
1355                 xi->ChangeMassHypothesis(lV0quality , 3312);    // Back to default hyp.
1356         }// end if negative bachelor
1357
1358
1359         if( bachTrackXi->Charge() >  0 ){
1360                 lV0quality = 0.;
1361                 xi->ChangeMassHypothesis(lV0quality , -3312);
1362                         // Calculate the effective mass of the Xi+ candidate.
1363                         // pdg code -3312 = Xi+
1364                 lInvMassXiPlus = xi->GetEffMassXi();
1365
1366                 lV0quality = 0.;
1367                 xi->ChangeMassHypothesis(lV0quality , -3334);
1368                         // Calculate the effective mass of the Xi+ candidate.
1369                         // pdg code -3334  = Omega+
1370                 lInvMassOmegaPlus = xi->GetEffMassXi();
1371
1372                 lV0quality = 0.;
1373                 xi->ChangeMassHypothesis(lV0quality , -3312);   // Back to "default" hyp.
1374         }// end if positive bachelor
1375
1376
1377
1378                 // - II.Step 5 : PID on the daughter tracks
1379                 //-------------
1380
1381         // A - Combined PID
1382         // Reasonable guess for the priors for the cascade track sample (e-, mu, pi, K, p)
1383         Double_t lPriorsGuessXi[5]    = {0, 0, 2, 0, 1};
1384         Double_t lPriorsGuessOmega[5] = {0, 0, 1, 1, 1};
1385
1386         // Combined VO-positive-daughter PID
1387         AliPID pPidXi;          pPidXi.SetPriors(    lPriorsGuessXi    );
1388         AliPID pPidOmega;       pPidOmega.SetPriors( lPriorsGuessOmega );
1389
1390         if( pTrackXi->IsOn(AliESDtrack::kESDpid) ){  // Combined PID exists
1391                 Double_t r[10] = {0.}; pTrackXi->GetESDpid(r);
1392                 pPidXi.SetProbabilities(r);
1393                 pPidOmega.SetProbabilities(r);
1394
1395                 // Check if the V0 positive track is a proton (case for Xi-)
1396                 Double_t pproton = pPidXi.GetProbability(AliPID::kProton);
1397                 if (pproton > pPidXi.GetProbability(AliPID::kElectron) &&
1398                     pproton > pPidXi.GetProbability(AliPID::kMuon)     &&
1399                     pproton > pPidXi.GetProbability(AliPID::kPion)     &&
1400                     pproton > pPidXi.GetProbability(AliPID::kKaon)     )     lIsPosInXiProton = kTRUE;
1401
1402                 // Check if the V0 positive track is a pi+ (case for Xi+)
1403                 Double_t ppion = pPidXi.GetProbability(AliPID::kPion);
1404                 if (ppion > pPidXi.GetProbability(AliPID::kElectron) &&
1405                     ppion > pPidXi.GetProbability(AliPID::kMuon)     &&
1406                     ppion > pPidXi.GetProbability(AliPID::kKaon)     &&
1407                     ppion > pPidXi.GetProbability(AliPID::kProton)   )     lIsPosInXiPion = kTRUE;
1408
1409
1410                 // Check if the V0 positive track is a proton (case for Omega-)
1411                 pproton = 0.;
1412                     pproton = pPidOmega.GetProbability(AliPID::kProton);
1413                 if (pproton > pPidOmega.GetProbability(AliPID::kElectron) &&
1414                     pproton > pPidOmega.GetProbability(AliPID::kMuon)     &&
1415                     pproton > pPidOmega.GetProbability(AliPID::kPion)     &&
1416                     pproton > pPidOmega.GetProbability(AliPID::kKaon)     )  lIsPosInOmegaProton = kTRUE;
1417
1418                 // Check if the V0 positive track is a pi+ (case for Omega+)
1419                 ppion = 0.;
1420                     ppion = pPidOmega.GetProbability(AliPID::kPion);
1421                 if (ppion > pPidOmega.GetProbability(AliPID::kElectron) &&
1422                     ppion > pPidOmega.GetProbability(AliPID::kMuon)     &&
1423                     ppion > pPidOmega.GetProbability(AliPID::kKaon)     &&
1424                     ppion > pPidOmega.GetProbability(AliPID::kProton)   )    lIsPosInOmegaPion = kTRUE;
1425
1426         }// end if V0 positive track with existing combined PID
1427
1428
1429         // Combined VO-negative-daughter PID
1430         AliPID nPidXi;          nPidXi.SetPriors(    lPriorsGuessXi    );
1431         AliPID nPidOmega;       nPidOmega.SetPriors( lPriorsGuessOmega );
1432
1433         if( nTrackXi->IsOn(AliESDtrack::kESDpid) ){  // Combined PID exists
1434                 Double_t r[10] = {0.}; nTrackXi->GetESDpid(r);
1435                 nPidXi.SetProbabilities(r);
1436                 nPidOmega.SetProbabilities(r);
1437
1438                 // Check if the V0 negative track is a pi- (case for Xi-)
1439                 Double_t ppion = nPidXi.GetProbability(AliPID::kPion);
1440                 if (ppion > nPidXi.GetProbability(AliPID::kElectron) &&
1441                     ppion > nPidXi.GetProbability(AliPID::kMuon)     &&
1442                     ppion > nPidXi.GetProbability(AliPID::kKaon)     &&
1443                     ppion > nPidXi.GetProbability(AliPID::kProton)   )     lIsNegInXiPion = kTRUE;
1444
1445                 // Check if the V0 negative track is an anti-proton (case for Xi+)
1446                 Double_t pproton = nPidXi.GetProbability(AliPID::kProton);
1447                 if (pproton > nPidXi.GetProbability(AliPID::kElectron) &&
1448                     pproton > nPidXi.GetProbability(AliPID::kMuon)     &&
1449                     pproton > nPidXi.GetProbability(AliPID::kPion)     &&
1450                     pproton > nPidXi.GetProbability(AliPID::kKaon)     )     lIsNegInXiProton = kTRUE;
1451
1452                 // Check if the V0 negative track is a pi- (case for Omega-)
1453                 ppion = 0.;
1454                     ppion = nPidOmega.GetProbability(AliPID::kPion);
1455                 if (ppion > nPidOmega.GetProbability(AliPID::kElectron) &&
1456                     ppion > nPidOmega.GetProbability(AliPID::kMuon)     &&
1457                     ppion > nPidOmega.GetProbability(AliPID::kKaon)     &&
1458                     ppion > nPidOmega.GetProbability(AliPID::kProton)   )    lIsNegInOmegaPion = kTRUE;
1459
1460                 // Check if the V0 negative track is an anti-proton (case for Omega+)
1461                 pproton = 0.;
1462                          pproton = nPidOmega.GetProbability(AliPID::kProton);
1463                 if (pproton > nPidOmega.GetProbability(AliPID::kElectron) &&
1464                     pproton > nPidOmega.GetProbability(AliPID::kMuon)     &&
1465                     pproton > nPidOmega.GetProbability(AliPID::kPion)     &&
1466                     pproton > nPidOmega.GetProbability(AliPID::kKaon)     )  lIsNegInOmegaProton = kTRUE;
1467
1468         }// end if V0 negative track with existing combined PID
1469
1470
1471         // Combined bachelor PID
1472         AliPID bachPidXi;       bachPidXi.SetPriors(    lPriorsGuessXi    );
1473         AliPID bachPidOmega;    bachPidOmega.SetPriors( lPriorsGuessOmega );
1474
1475         if( bachTrackXi->IsOn(AliESDtrack::kESDpid) ){  // Combined PID exists
1476                 Double_t r[10] = {0.}; bachTrackXi->GetESDpid(r);
1477                 bachPidXi.SetProbabilities(r);
1478                 bachPidOmega.SetProbabilities(r);
1479                 // Check if the bachelor track is a pion
1480                 Double_t ppion = bachPidXi.GetProbability(AliPID::kPion);
1481                 if (ppion > bachPidXi.GetProbability(AliPID::kElectron) &&
1482                     ppion > bachPidXi.GetProbability(AliPID::kMuon)     &&
1483                     ppion > bachPidXi.GetProbability(AliPID::kKaon)     &&
1484                     ppion > bachPidXi.GetProbability(AliPID::kProton)   )     lIsBachelorPion = kTRUE;
1485                 // Check if the bachelor track is a kaon
1486                 Double_t pkaon = bachPidOmega.GetProbability(AliPID::kKaon);
1487                 if (pkaon > bachPidOmega.GetProbability(AliPID::kElectron) &&
1488                     pkaon > bachPidOmega.GetProbability(AliPID::kMuon)     &&
1489                     pkaon > bachPidOmega.GetProbability(AliPID::kPion)     &&
1490                     pkaon > bachPidOmega.GetProbability(AliPID::kProton)   )  lIsBachelorKaon = kTRUE;
1491         }// end if bachelor track with existing combined PID
1492
1493
1494         // B - TPC PID : 3-sigma bands on Bethe-Bloch curve
1495         // Bachelor
1496         if (TMath::Abs(fTpcPidManager->GetNumberOfSigmas(bachTrackXi,AliPID::kKaon)) < 3) lIsBachelorKaonForTPC = kTRUE;
1497         if (TMath::Abs(fTpcPidManager->GetNumberOfSigmas(bachTrackXi,AliPID::kPion)) < 3) lIsBachelorPionForTPC = kTRUE;
1498
1499         // Negative V0 daughter
1500         if (TMath::Abs(fTpcPidManager->GetNumberOfSigmas(nTrackXi,AliPID::kPion   )) < 3) lIsNegPionForTPC   = kTRUE;
1501         if (TMath::Abs(fTpcPidManager->GetNumberOfSigmas(nTrackXi,AliPID::kProton )) < 3) lIsNegProtonForTPC = kTRUE;
1502
1503         // Positive V0 daughter
1504         if (TMath::Abs(fTpcPidManager->GetNumberOfSigmas(pTrackXi,AliPID::kPion   )) < 3) lIsPosPionForTPC   = kTRUE;
1505         if (TMath::Abs(fTpcPidManager->GetNumberOfSigmas(pTrackXi,AliPID::kProton )) < 3) lIsPosProtonForTPC = kTRUE;
1506
1507
1508
1509                 // - II.Step 6 : extra info for QA (ESD)
1510                 // miscellaneous pieces of info that may help regarding data quality assessment.
1511                 //-------------
1512
1513         xi->GetPxPyPz( lXiMomX, lXiMomY, lXiMomZ );
1514                 lXiTransvMom    = TMath::Sqrt( lXiMomX*lXiMomX   + lXiMomY*lXiMomY );
1515                 lXiTotMom       = TMath::Sqrt( lXiMomX*lXiMomX   + lXiMomY*lXiMomY   + lXiMomZ*lXiMomZ );
1516
1517         xi->GetBPxPyPz(  lBachMomX,  lBachMomY,  lBachMomZ );
1518                 lBachTransvMom  = TMath::Sqrt( lBachMomX*lBachMomX   + lBachMomY*lBachMomY );
1519                 lBachTotMom     = TMath::Sqrt( lBachMomX*lBachMomX   + lBachMomY*lBachMomY  +  lBachMomZ*lBachMomZ  );
1520
1521         lChargeXi = xi->Charge();
1522
1523         lV0toXiCosineOfPointingAngle = xi->GetV0CosineOfPointingAngle( lPosXi[0], lPosXi[1], lPosXi[2] );
1524
1525         lRapXi    = xi->RapXi();
1526         lRapOmega = xi->RapOmega();
1527         lEta      = xi->Eta();
1528         lTheta    = xi->Theta() *180.0/TMath::Pi();
1529         lPhi      = xi->Phi()   *180.0/TMath::Pi();
1530         lAlphaXi  = xi->AlphaXi();
1531         lPtArmXi  = xi->PtArmXi();
1532
1533
1534
1535                 // II.Step 7 - Complementary info for monitoring the cascade cut variables
1536
1537         const AliMultiplicity *lAliMult = lESDevent->GetMultiplicity();
1538         lSPDTrackletsMultiplicity = lAliMult->GetNumberOfTracklets();
1539         lBachTPCClusters          = bachTrackXi->GetTPCNcls(); // for the Bachelor at least, (we will see later for the baryon, for the meson)
1540
1541
1542                 // II.Step 8 - Azimuthal correlation study
1543                 //-------------
1544
1545         lTVect3MomXi.SetXYZ( lXiMomX, lXiMomY, lXiMomZ );
1546         lArrTrackID[0] = pTrackXi   ->GetID();
1547         lArrTrackID[1] = nTrackXi   ->GetID();
1548         lArrTrackID[2] = bachTrackXi->GetID();
1549
1550
1551   }// end of ESD treatment
1552
1553
1554   if(fAnalysisType == "AOD"){
1555
1556         // -------------------------------------
1557         // II.AOD - Calcultaion Part dedicated to Xi vertices (ESD)
1558
1559         const AliAODcascade *xi = lAODevent->GetCascade(iXi);
1560         if (!xi) continue;
1561
1562         // Just to know which file is currently open : locate the file containing Xi
1563         // cout << "Name of the file containing Xi candidate(s) :" <<  fesdH->GetTree()->GetCurrentFile()->GetName() << endl;
1564
1565
1566                 // - II.Step 1 : Characteristics of the event : prim. Vtx + magnetic field (AOD)
1567                 //-------------
1568
1569         lTrkgPrimaryVtxPos[0]   = -100.0;
1570         lTrkgPrimaryVtxPos[1]   = -100.0;
1571         lTrkgPrimaryVtxPos[2]   = -100.0;
1572         lTrkgPrimaryVtxRadius3D = -500. ;
1573         // We don't have the different prim. vertex at the AOD level -> nothing to do.
1574
1575         const AliAODVertex *lPrimaryBestVtx = lAODevent->GetPrimaryVertex();
1576         // get the best primary vertex available for the event
1577         // We may keep the one which is the best one available = GetVertex(0)
1578         // Pb with pile-up to expect
1579         // This one will be used for next calculations (DCA essentially)
1580
1581         lPrimaryBestVtx->GetXYZ( lBestPrimaryVtxPos );
1582         lBestPrimaryVtxRadius3D = TMath::Sqrt(  lBestPrimaryVtxPos[0] * lBestPrimaryVtxPos[0] +
1583                                                 lBestPrimaryVtxPos[1] * lBestPrimaryVtxPos[1] +
1584                                                 lBestPrimaryVtxPos[2] * lBestPrimaryVtxPos[2] );
1585
1586
1587                 // - II.Step 2 : Assigning the necessary variables for specific AliAODcascade data members (AOD)
1588                 //-------------
1589
1590         lEffMassXi                      = xi->MassXi(); // default working hypothesis : cascade = Xi- decay
1591         lChi2Xi                         = xi->Chi2Xi();
1592         lDcaXiDaughters                 = xi->DcaXiDaughters();
1593         lXiCosineOfPointingAngle        = xi->CosPointingAngleXi( lBestPrimaryVtxPos[0],
1594                                                                   lBestPrimaryVtxPos[1],
1595                                                                   lBestPrimaryVtxPos[2] );
1596                                         // Take care :
1597                                         // the best available vertex should be used (like in AliCascadeVertexer)
1598
1599                 lPosXi[0] = xi->DecayVertexXiX();
1600                 lPosXi[1] = xi->DecayVertexXiY();
1601                 lPosXi[2] = xi->DecayVertexXiZ();
1602         lXiRadius = TMath::Sqrt( lPosXi[0]*lPosXi[0]  +  lPosXi[1]*lPosXi[1] );
1603
1604
1605                 // - II.Step 3 : around the tracks : Bach + V0 (AOD)
1606                 // ~ Necessary variables for AODcascade data members coming from the AODv0 part (inheritance)
1607                 //-------------
1608
1609         lChargeXi                       = xi->ChargeXi();
1610
1611         if( lChargeXi < 0)
1612           lInvMassLambdaAsCascDghter    = xi->MassLambda();
1613         else
1614           lInvMassLambdaAsCascDghter    = xi->MassAntiLambda();
1615
1616         lDcaV0DaughtersXi               = xi->DcaV0Daughters();
1617         lV0Chi2Xi                       = xi->Chi2V0();
1618
1619         lV0CosineOfPointingAngleXi      = xi->CosPointingAngle( lBestPrimaryVtxPos );
1620         lDcaV0ToPrimVertexXi            = xi->DcaV0ToPrimVertex();
1621
1622         lDcaBachToPrimVertexXi          = xi->DcaBachToPrimVertex();
1623
1624
1625                 lPosV0Xi[0] = xi->DecayVertexV0X();
1626                 lPosV0Xi[1] = xi->DecayVertexV0Y();
1627                 lPosV0Xi[2] = xi->DecayVertexV0Z();
1628         lV0RadiusXi     = TMath::Sqrt( lPosV0Xi[0]*lPosV0Xi[0]  +  lPosV0Xi[1]*lPosV0Xi[1] );
1629
1630         lDcaPosToPrimVertexXi           = xi->DcaPosToPrimVertex();
1631         lDcaNegToPrimVertexXi           = xi->DcaNegToPrimVertex();
1632
1633
1634                 // - II.Step 4 : around effective masses (AOD)
1635                 // ~ change mass hypotheses to cover all the possibilities :  Xi-/+, Omega -/+
1636                 //-------------
1637
1638         if( lChargeXi < 0 )             lInvMassXiMinus         = xi->MassXi();
1639         if( lChargeXi > 0 )             lInvMassXiPlus          = xi->MassXi();
1640         if( lChargeXi < 0 )             lInvMassOmegaMinus      = xi->MassOmega();
1641         if( lChargeXi > 0 )             lInvMassOmegaPlus       = xi->MassOmega();
1642
1643
1644                 // - II.Step 5 : PID on the daughters (To be developed ...)
1645                 //-------------
1646
1647         // Combined PID
1648
1649         /*
1650         // Reasonable guess for the priors for the cascade track sample
1651         Double_t lPriorsGuessXi[5]    = {0.0, 0.0, 2, 0, 1};
1652         Double_t lPriorsGuessOmega[5] = {0.0, 0.0, 1, 1, 1};
1653         AliPID bachPidXi;       bachPidXi.SetPriors(    lPriorsGuessXi    );
1654         AliPID bachPidOmega;    bachPidOmega.SetPriors( lPriorsGuessOmega );
1655
1656         const AliAODTrack *bachTrackXi = lAODevent->GetTrack( xi->GetBachID() ); // FIXME : GetBachID not implemented ?
1657
1658         if( bachTrackXi->IsOn(AliESDtrack::kESDpid) ){  // Combined PID exists, the AOD flags = a copy of the ESD ones
1659                 Double_t r[10]; bachTrackXi->GetPID(r);
1660                 bachPidXi.SetProbabilities(r);
1661                 bachPidOmega.SetProbabilities(r);
1662                 // Check if the bachelor track is a pion
1663                 Double_t ppion = bachPidXi.GetProbability(AliPID::kPion);
1664                 if (ppion > bachPidXi.GetProbability(AliPID::kElectron) &&
1665                     ppion > bachPidXi.GetProbability(AliPID::kMuon)     &&
1666                     ppion > bachPidXi.GetProbability(AliPID::kKaon)     &&
1667                     ppion > bachPidXi.GetProbability(AliPID::kProton)   )     lIsBachelorPion = kTRUE;
1668                 // Check if the bachelor track is a kaon
1669                 Double_t pkaon = bachPidOmega.GetProbability(AliPID::kKaon);
1670                 if (pkaon > bachPidOmega.GetProbability(AliPID::kElectron) &&
1671                     pkaon > bachPidOmega.GetProbability(AliPID::kMuon)     &&
1672                     pkaon > bachPidOmega.GetProbability(AliPID::kPion)     &&
1673                     pkaon > bachPidOmega.GetProbability(AliPID::kProton)   )  lIsBachelorKaon = kTRUE;
1674
1675         }// end if bachelor track with existing combined PID
1676         */
1677
1678         // TPC PID
1679
1680                 // - II.Step 6 : extra info for QA (AOD)
1681                 // miscellaneous pieces onf info that may help regarding data quality assessment.
1682                 //-------------
1683
1684                 lXiMomX = xi->MomXiX();
1685                 lXiMomY = xi->MomXiY();
1686                 lXiMomZ = xi->MomXiZ();
1687         lXiTransvMom    = TMath::Sqrt( lXiMomX*lXiMomX   + lXiMomY*lXiMomY );
1688         lXiTotMom       = TMath::Sqrt( lXiMomX*lXiMomX   + lXiMomY*lXiMomY   + lXiMomZ*lXiMomZ );
1689
1690                 lBachMomX = xi->MomBachX();
1691                 lBachMomY = xi->MomBachY();
1692                 lBachMomZ = xi->MomBachZ();
1693         lBachTransvMom  = TMath::Sqrt( lBachMomX*lBachMomX   + lBachMomY*lBachMomY );
1694         lBachTotMom     = TMath::Sqrt( lBachMomX*lBachMomX   + lBachMomY*lBachMomY  +  lBachMomZ*lBachMomZ  );
1695
1696
1697         lV0toXiCosineOfPointingAngle = xi->CosPointingAngle( xi->GetDecayVertexXi() );
1698
1699         lRapXi    = xi->RapXi();
1700         lRapOmega = xi->RapOmega();
1701         lEta      = xi->Eta();                          // Will not work ! need a method Pz(), Py() Px()
1702         lTheta    = xi->Theta() *180.0/TMath::Pi();     // in AODcascade.
1703         lPhi      = xi->Phi()   *180.0/TMath::Pi();     // Here, we will get eta, theta, phi for the V0 ...
1704         lAlphaXi  = xi->AlphaXi();
1705         lPtArmXi  = xi->PtArmXi();
1706
1707                 // II.Step 7 - Complementary info for monitoring the cascade cut variables
1708         //FIXME
1709
1710                 // II.Step 8 - Azimuthal correlation study
1711                 //-------------
1712
1713         lTVect3MomXi.SetXYZ( lXiMomX, lXiMomY, lXiMomZ );
1714
1715         AliAODTrack *pTrackXi    = dynamic_cast<AliAODTrack*>( xi->GetDaughter(0) );
1716         AliAODTrack *nTrackXi    = dynamic_cast<AliAODTrack*>( xi->GetDaughter(1) );
1717         AliAODTrack *bachTrackXi = dynamic_cast<AliAODTrack*>( xi->GetDecayVertexXi()->GetDaughter(0) );
1718                 if (!pTrackXi || !nTrackXi || !bachTrackXi ) {
1719                         Printf("ERROR: Could not retrieve one of the 3 AOD daughter tracks of the cascade ...");
1720                         continue;
1721                 }
1722
1723         lArrTrackID[0] = pTrackXi   ->GetID();
1724         lArrTrackID[1] = nTrackXi   ->GetID();
1725         lArrTrackID[2] = bachTrackXi->GetID();
1726
1727   }// end of AOD treatment
1728
1729
1730     // -------------------------------------
1731     // II.Fill - Filling the TH1,2,3Fs
1732
1733
1734         // - II.Fill.Step 1      : primary vertex
1735
1736         fHistVtxStatus                  ->Fill( lStatusTrackingPrimVtx   );  // 1 if tracking vtx = ok
1737
1738         if( lStatusTrackingPrimVtx ){
1739                 fHistPosTrkgPrimaryVtxX ->Fill( lTrkgPrimaryVtxPos[0]    );
1740                 fHistPosTrkgPrimaryVtxY ->Fill( lTrkgPrimaryVtxPos[1]    );
1741                 fHistPosTrkgPrimaryVtxZ ->Fill( lTrkgPrimaryVtxPos[2]    );
1742                 fHistTrkgPrimaryVtxRadius->Fill( lTrkgPrimaryVtxRadius3D );
1743         }
1744
1745         fHistPosBestPrimaryVtxX         ->Fill( lBestPrimaryVtxPos[0]    );
1746         fHistPosBestPrimaryVtxY         ->Fill( lBestPrimaryVtxPos[1]    );
1747         fHistPosBestPrimaryVtxZ         ->Fill( lBestPrimaryVtxPos[2]    );
1748         fHistBestPrimaryVtxRadius       ->Fill( lBestPrimaryVtxRadius3D  );
1749
1750         f2dHistTrkgPrimVtxVsBestPrimVtx->Fill( lTrkgPrimaryVtxRadius3D, lBestPrimaryVtxRadius3D );
1751
1752
1753         // II.Fill.Step 2
1754         fHistEffMassXi                  ->Fill( lEffMassXi               );
1755         fHistChi2Xi                     ->Fill( lChi2Xi                  );     // Flag CascadeVtxer: Cut Variable a
1756         fHistDcaXiDaughters             ->Fill( lDcaXiDaughters          );     // Flag CascadeVtxer: Cut Variable e
1757         fHistDcaBachToPrimVertex        ->Fill( lDcaBachToPrimVertexXi   );     // Flag CascadeVtxer: Cut Variable d
1758         fHistXiCosineOfPointingAngle    ->Fill( lXiCosineOfPointingAngle );     // Flag CascadeVtxer: Cut Variable f
1759         fHistXiRadius                   ->Fill( lXiRadius                );     // Flag CascadeVtxer: Cut Variable g+h
1760
1761
1762         // II.Fill.Step 3
1763         fHistMassLambdaAsCascDghter     ->Fill( lInvMassLambdaAsCascDghter );   // Flag CascadeVtxer: Cut Variable c
1764         fHistV0Chi2Xi                   ->Fill( lV0Chi2Xi                  );
1765         fHistDcaV0DaughtersXi           ->Fill( lDcaV0DaughtersXi          );
1766         fHistV0CosineOfPointingAngleXi  ->Fill( lV0CosineOfPointingAngleXi );
1767         fHistV0RadiusXi                 ->Fill( lV0RadiusXi                );
1768
1769         fHistDcaV0ToPrimVertexXi        ->Fill( lDcaV0ToPrimVertexXi       );   // Flag CascadeVtxer: Cut Variable b
1770         fHistDcaPosToPrimVertexXi       ->Fill( lDcaPosToPrimVertexXi      );
1771         fHistDcaNegToPrimVertexXi       ->Fill( lDcaNegToPrimVertexXi      );
1772
1773
1774         // II.Fill.Step 4 : extra QA info
1775         fHistXiTransvMom        ->Fill( lXiTransvMom   );
1776         fHistXiTotMom           ->Fill( lXiTotMom      );
1777
1778         fHistBachTransvMom      ->Fill( lBachTransvMom );
1779         fHistBachTotMom         ->Fill( lBachTotMom    );
1780
1781         fHistChargeXi           ->Fill( lChargeXi      );
1782         fHistV0toXiCosineOfPointingAngle->Fill( lV0toXiCosineOfPointingAngle );
1783
1784         fHistRapXi              ->Fill( lRapXi               );
1785         fHistRapOmega           ->Fill( lRapOmega            );
1786         fHistEta                ->Fill( lEta                 );
1787         fHistTheta              ->Fill( lTheta               );
1788         fHistPhi                ->Fill( lPhi                 );
1789
1790         f2dHistArmenteros       ->Fill( lAlphaXi, lPtArmXi   );
1791
1792
1793         // II.Fill.Step 5 : inv mass plots 1D
1794         if( lChargeXi < 0 ){
1795                                         fHistMassXiMinus               ->Fill( lInvMassXiMinus    );
1796                                         fHistMassOmegaMinus            ->Fill( lInvMassOmegaMinus );
1797                 if(lIsBachelorPion)     fHistMassWithCombPIDXiMinus    ->Fill( lInvMassXiMinus    );
1798                 if(lIsBachelorKaon)     fHistMassWithCombPIDOmegaMinus ->Fill( lInvMassOmegaMinus );
1799         }
1800
1801         if( lChargeXi > 0 ){
1802                                         fHistMassXiPlus                ->Fill( lInvMassXiPlus     );
1803                                         fHistMassOmegaPlus             ->Fill( lInvMassOmegaPlus  );
1804                 if(lIsBachelorPion)     fHistMassWithCombPIDXiPlus     ->Fill( lInvMassXiPlus     );
1805                 if(lIsBachelorKaon)     fHistMassWithCombPIDOmegaPlus  ->Fill( lInvMassOmegaPlus  );
1806         }
1807
1808
1809         // II.Fill.Step 6 : inv mass plots 2D, 3D
1810         if( lChargeXi < 0 ) {
1811                 f2dHistEffMassLambdaVsEffMassXiMinus->Fill( lInvMassLambdaAsCascDghter, lInvMassXiMinus );
1812                 f2dHistEffMassXiVsEffMassOmegaMinus ->Fill( lInvMassXiMinus, lInvMassOmegaMinus );
1813                 f2dHistXiRadiusVsEffMassXiMinus     ->Fill( lXiRadius, lInvMassXiMinus );
1814                 f2dHistXiRadiusVsEffMassOmegaMinus  ->Fill( lXiRadius, lInvMassOmegaMinus );
1815                 f3dHistXiPtVsEffMassVsYXiMinus      ->Fill( lXiTransvMom, lInvMassXiMinus,    lRapXi    );
1816                 f3dHistXiPtVsEffMassVsYOmegaMinus   ->Fill( lXiTransvMom, lInvMassOmegaMinus, lRapOmega );
1817                 if(lIsPosInXiProton)                            f3dHistXiPtVsEffMassVsYWithCombPIDXiMinus     ->Fill( lXiTransvMom, lInvMassXiMinus,    lRapXi    );
1818                 if(lIsBachelorKaon)                             f3dHistXiPtVsEffMassVsYWithCombPIDOmegaMinus  ->Fill( lXiTransvMom, lInvMassOmegaMinus, lRapOmega );
1819                 if(lIsBachelorKaon && lIsPosInOmegaProton)      f3dHistXiPtVsEffMassVsYWith2CombPIDOmegaMinus ->Fill( lXiTransvMom, lInvMassOmegaMinus, lRapOmega );
1820                 if(lIsBachelorKaonForTPC)                       f3dHistXiPtVsEffMassVsYWithTpcPIDOmegaMinus   ->Fill( lXiTransvMom, lInvMassOmegaMinus, lRapOmega );
1821         }
1822         else{
1823                 f2dHistEffMassLambdaVsEffMassXiPlus ->Fill( lInvMassLambdaAsCascDghter, lInvMassXiPlus );
1824                 f2dHistEffMassXiVsEffMassOmegaPlus  ->Fill( lInvMassXiPlus, lInvMassOmegaPlus );
1825                 f2dHistXiRadiusVsEffMassXiPlus      ->Fill( lXiRadius, lInvMassXiPlus);
1826                 f2dHistXiRadiusVsEffMassOmegaPlus   ->Fill( lXiRadius, lInvMassOmegaPlus );
1827                 f3dHistXiPtVsEffMassVsYXiPlus       ->Fill( lXiTransvMom, lInvMassXiPlus,    lRapXi    );
1828                 f3dHistXiPtVsEffMassVsYOmegaPlus    ->Fill( lXiTransvMom, lInvMassOmegaPlus, lRapOmega );
1829                 if(lIsNegInXiProton)                            f3dHistXiPtVsEffMassVsYWithCombPIDXiPlus     ->Fill( lXiTransvMom, lInvMassXiPlus,    lRapXi    );
1830                 if(lIsBachelorKaon )                            f3dHistXiPtVsEffMassVsYWithCombPIDOmegaPlus  ->Fill( lXiTransvMom, lInvMassOmegaPlus, lRapOmega );
1831                 if(lIsBachelorKaon && lIsNegInOmegaProton)      f3dHistXiPtVsEffMassVsYWith2CombPIDOmegaPlus ->Fill( lXiTransvMom, lInvMassOmegaPlus, lRapOmega );
1832         }
1833
1834         // - Filling the AliCFContainers related to PID
1835
1836         Double_t lContainerPIDVars[4] = {0.0};
1837
1838         // Xi Minus
1839         if( lChargeXi < 0 ) {
1840                 lContainerPIDVars[0] = lXiTransvMom       ;
1841                 lContainerPIDVars[1] = lInvMassXiMinus    ;
1842                 lContainerPIDVars[2] = lRapXi             ;
1843                 lContainerPIDVars[3] = nTrackMultiplicity ;
1844
1845                 // No PID
1846                         fCFContCascadePIDXiMinus->Fill(lContainerPIDVars, 0); // No PID
1847                 // TPC PID
1848                 if( lIsBachelorPionForTPC  )
1849                         fCFContCascadePIDXiMinus->Fill(lContainerPIDVars, 1); // TPC PID / 3-#sigma cut on Bachelor track
1850
1851                 if( lIsBachelorPionForTPC &&
1852                     lIsPosProtonForTPC     )
1853                         fCFContCascadePIDXiMinus->Fill(lContainerPIDVars, 2); // TPC PID / 3-#sigma cut on Bachelor+Baryon tracks
1854
1855                 if( lIsBachelorPionForTPC &&
1856                     lIsPosProtonForTPC    &&
1857                     lIsNegPionForTPC       )
1858                         fCFContCascadePIDXiMinus->Fill(lContainerPIDVars, 3); // TPC PID / 3-#sigma cut on Bachelor+Baryon+Meson tracks
1859
1860                 // Combined PID
1861                 if( lIsBachelorPion        )
1862                         fCFContCascadePIDXiMinus->Fill(lContainerPIDVars, 4); // Comb. PID / Bachelor
1863
1864                 if( lIsBachelorPion       &&
1865                     lIsPosInXiProton    )
1866                         fCFContCascadePIDXiMinus->Fill(lContainerPIDVars, 5); // Comb. PID / Bachelor+Baryon
1867
1868                 if(lIsBachelorPion     &&
1869                    lIsPosInXiProton &&
1870                    lIsNegInXiPion    )
1871                         fCFContCascadePIDXiMinus->Fill(lContainerPIDVars, 6); // Comb. PID / Bachelor+Baryon+Meson
1872         }
1873
1874         lContainerPIDVars[0] = 0.; lContainerPIDVars[1] = 0.; lContainerPIDVars[2] = 0.; lContainerPIDVars[3] = 0.;
1875
1876         // Xi Plus
1877         if( lChargeXi > 0 ) {
1878                 lContainerPIDVars[0] = lXiTransvMom       ;
1879                 lContainerPIDVars[1] = lInvMassXiPlus     ;
1880                 lContainerPIDVars[2] = lRapXi             ;
1881                 lContainerPIDVars[3] = nTrackMultiplicity ;
1882
1883                 // No PID
1884                         fCFContCascadePIDXiPlus->Fill(lContainerPIDVars, 0); // No PID
1885                 // TPC PID
1886                 if( lIsBachelorPionForTPC  )
1887                         fCFContCascadePIDXiPlus->Fill(lContainerPIDVars, 1); // TPC PID / 3-#sigma cut on Bachelor track
1888
1889                 if( lIsBachelorPionForTPC &&
1890                     lIsNegProtonForTPC     )
1891                         fCFContCascadePIDXiPlus->Fill(lContainerPIDVars, 2); // TPC PID / 3-#sigma cut on Bachelor+Baryon tracks
1892
1893                 if( lIsBachelorPionForTPC &&
1894                     lIsNegProtonForTPC    &&
1895                     lIsPosPionForTPC       )
1896                         fCFContCascadePIDXiPlus->Fill(lContainerPIDVars, 3); // TPC PID / 3-#sigma cut on Bachelor+Baryon+Meson tracks
1897
1898                 // Combined PID
1899                 if( lIsBachelorPion        )
1900                         fCFContCascadePIDXiPlus->Fill(lContainerPIDVars, 4); // Comb. PID / Bachelor
1901
1902                 if( lIsBachelorPion       &&
1903                     lIsNegInXiProton    )
1904                         fCFContCascadePIDXiPlus->Fill(lContainerPIDVars, 5); // Comb. PID / Bachelor+Baryon
1905
1906                 if(lIsBachelorPion     &&
1907                    lIsNegInXiProton &&
1908                    lIsPosInXiPion    )
1909                         fCFContCascadePIDXiPlus->Fill(lContainerPIDVars, 6); // Comb. PID / Bachelor+Baryon+Meson
1910         }
1911
1912         lContainerPIDVars[0] = 0.; lContainerPIDVars[1] = 0.; lContainerPIDVars[2] = 0.; lContainerPIDVars[3] = 0.;
1913
1914         // Omega Minus
1915         if( lChargeXi < 0 ) {
1916                 lContainerPIDVars[0] = lXiTransvMom       ;
1917                 lContainerPIDVars[1] = lInvMassOmegaMinus ;
1918                 lContainerPIDVars[2] = lRapOmega          ;
1919                 lContainerPIDVars[3] = nTrackMultiplicity ;
1920
1921                 // No PID
1922                         fCFContCascadePIDOmegaMinus->Fill(lContainerPIDVars, 0); // No PID
1923                 // TPC PID
1924                 if( lIsBachelorKaonForTPC  )
1925                         fCFContCascadePIDOmegaMinus->Fill(lContainerPIDVars, 1); // TPC PID / 3-#sigma cut on Bachelor track
1926
1927                 if( lIsBachelorKaonForTPC &&
1928                     lIsPosProtonForTPC     )
1929                         fCFContCascadePIDOmegaMinus->Fill(lContainerPIDVars, 2); // TPC PID / 3-#sigma cut on Bachelor+Baryon tracks
1930
1931                 if( lIsBachelorKaonForTPC &&
1932                     lIsPosProtonForTPC    &&
1933                     lIsNegPionForTPC       )
1934                         fCFContCascadePIDOmegaMinus->Fill(lContainerPIDVars, 3); // TPC PID / 3-#sigma cut on Bachelor+Baryon+Meson tracks
1935
1936                 // Combined PID
1937                 if( lIsBachelorKaon        )
1938                         fCFContCascadePIDOmegaMinus->Fill(lContainerPIDVars, 4); // Comb. PID / Bachelor
1939
1940                 if( lIsBachelorKaon       &&
1941                     lIsPosInOmegaProton    )
1942                         fCFContCascadePIDOmegaMinus->Fill(lContainerPIDVars, 5); // Comb. PID / Bachelor+Baryon
1943
1944                 if(lIsBachelorKaon     &&
1945                    lIsPosInOmegaProton &&
1946                    lIsNegInOmegaPion    )
1947                         fCFContCascadePIDOmegaMinus->Fill(lContainerPIDVars, 6); // Comb. PID / Bachelor+Baryon+Meson
1948         }
1949
1950         lContainerPIDVars[0] = 0.; lContainerPIDVars[1] = 0.; lContainerPIDVars[2] = 0.; lContainerPIDVars[3] = 0.;
1951
1952         // Omega Plus
1953         if( lChargeXi > 0 ) {
1954                 lContainerPIDVars[0] = lXiTransvMom       ;
1955                 lContainerPIDVars[1] = lInvMassOmegaPlus ;
1956                 lContainerPIDVars[2] = lRapOmega          ;
1957                 lContainerPIDVars[3] = nTrackMultiplicity ;
1958
1959                 // No PID
1960                         fCFContCascadePIDOmegaPlus->Fill(lContainerPIDVars, 0); // No PID
1961                 // TPC PID
1962                 if( lIsBachelorKaonForTPC  )
1963                         fCFContCascadePIDOmegaPlus->Fill(lContainerPIDVars, 1); // TPC PID / 3-#sigma cut on Bachelor track
1964
1965                 if( lIsBachelorKaonForTPC &&
1966                     lIsNegProtonForTPC     )
1967                         fCFContCascadePIDOmegaPlus->Fill(lContainerPIDVars, 2); // TPC PID / 3-#sigma cut on Bachelor+Baryon tracks
1968
1969                 if( lIsBachelorKaonForTPC &&
1970                     lIsNegProtonForTPC    &&
1971                     lIsPosPionForTPC       )
1972                         fCFContCascadePIDOmegaPlus->Fill(lContainerPIDVars, 3); // TPC PID / 3-#sigma cut on Bachelor+Baryon+Meson tracks
1973
1974                 // Combined PID
1975                 if( lIsBachelorKaon        )
1976                         fCFContCascadePIDOmegaPlus->Fill(lContainerPIDVars, 4); // Comb. PID / Bachelor
1977
1978                 if( lIsBachelorKaon       &&
1979                     lIsNegInOmegaProton    )
1980                         fCFContCascadePIDOmegaPlus->Fill(lContainerPIDVars, 5); // Comb. PID / Bachelor+Baryon
1981
1982                 if(lIsBachelorKaon     &&
1983                    lIsNegInOmegaProton &&
1984                    lIsPosInOmegaPion    )
1985                         fCFContCascadePIDOmegaPlus->Fill(lContainerPIDVars, 6); // Comb. PID / Bachelor+Baryon+Meson
1986         }
1987
1988
1989
1990
1991
1992         // II.Fill.Step 7 : filling the AliCFContainer (optimisation of topological selections)
1993         Double_t lContainerCutVars[18] = {0.0};
1994
1995         lContainerCutVars[0]  = lDcaXiDaughters;
1996         lContainerCutVars[1]  = lDcaBachToPrimVertexXi;
1997         lContainerCutVars[2]  = lXiCosineOfPointingAngle;
1998         lContainerCutVars[3]  = lXiRadius;
1999         lContainerCutVars[4]  = lInvMassLambdaAsCascDghter;
2000         lContainerCutVars[5]  = lDcaV0DaughtersXi;
2001         lContainerCutVars[6]  = lV0CosineOfPointingAngleXi;
2002         lContainerCutVars[7]  = lV0RadiusXi;
2003         lContainerCutVars[8]  = lDcaV0ToPrimVertexXi;
2004         lContainerCutVars[9]  = lDcaPosToPrimVertexXi;
2005         lContainerCutVars[10] = lDcaNegToPrimVertexXi;
2006
2007         if( lChargeXi < 0 ) {
2008                 lContainerCutVars[11] = lInvMassXiMinus;
2009                 lContainerCutVars[12] = lInvMassOmegaMinus;
2010         }
2011         else {
2012                 lContainerCutVars[11] = lInvMassXiPlus;
2013                 lContainerCutVars[12] = lInvMassOmegaPlus;
2014         }
2015
2016         lContainerCutVars[13] = lXiTransvMom;
2017         lContainerCutVars[14] = lBestPrimaryVtxPos[2];
2018         lContainerCutVars[15] = nTrackMultiplicity;
2019         lContainerCutVars[16] = lSPDTrackletsMultiplicity; // FIXME : SPDTrackletsMultiplicity is not available for AOD ... = -1
2020         lContainerCutVars[17] = lBachTPCClusters;          // FIXME : BachTPCClusters          is not available for AOD ... = -1
2021
2022         if( lChargeXi < 0 ) fCFContCascadeCuts->Fill(lContainerCutVars,0); // for negative cascades = Xi- and Omega-
2023         else                fCFContCascadeCuts->Fill(lContainerCutVars,1); // for negative cascades = Xi+ and Omega+
2024
2025
2026         // II.Fill.Step 8 :  angular correlations
2027
2028         if( lChargeXi < 0 ){
2029                 DoAngularCorrelation("Xi-",    lInvMassXiMinus,    lArrTrackID, lTVect3MomXi, lEta);
2030                 DoAngularCorrelation("Omega-", lInvMassOmegaMinus, lArrTrackID, lTVect3MomXi, lEta);
2031         }
2032         else{
2033                 DoAngularCorrelation("Xi+",    lInvMassXiPlus,    lArrTrackID, lTVect3MomXi, lEta);
2034                 DoAngularCorrelation("Omega+", lInvMassOmegaPlus, lArrTrackID, lTVect3MomXi, lEta);
2035         }
2036
2037
2038     }// end of the Cascade loop (ESD or AOD)
2039
2040
2041   // Post output data.
2042  PostData(1, fListHistCascade);
2043 }
2044
2045
2046 void AliAnalysisTaskCheckCascade::DoAngularCorrelation( const Char_t   *lCascType,
2047                                                               Double_t  lInvMassCascade,
2048                                                         const Int_t    *lArrTrackID,
2049                                                               TVector3 &lTVect3MomXi,
2050                                                               Double_t  lEtaXi ){
2051   // Perform the Delta(Phi)Delta(Eta) analysis
2052   // by properly filling the THnSparseF
2053
2054         TString lStrCascType( lCascType );
2055
2056         Double_t lCascPdgMass = 0.0;
2057         if( lStrCascType.Contains("Xi") )       lCascPdgMass = 1.3217;
2058         if( lStrCascType.Contains("Omega") )    lCascPdgMass = 1.6724;
2059
2060         if( lInvMassCascade > lCascPdgMass + 0.010) return;
2061         if( lInvMassCascade < lCascPdgMass - 0.010) return;
2062         // Check the Xi- candidate is within the proper mass window m0 +- 10 MeV
2063
2064
2065         // 1st loop: check there is no track with a higher pt ...
2066         // = The cascade is meant to be a leading particle : Pt(Casc) > any track in the event
2067         for(Int_t TrckIdx = 0; TrckIdx < (InputEvent())->GetNumberOfTracks() ; TrckIdx++ )
2068         {// Loop over all the tracks of the event
2069
2070                 AliVTrack *lCurrentTrck = dynamic_cast<AliVTrack*>( (InputEvent())->GetTrack( TrckIdx ) );
2071                         if (!lCurrentTrck ) {
2072                                 Printf("ERROR Correl. Study : Could not retrieve a track while looping over the event tracks ...");
2073                                 continue;
2074                         }
2075                 if(lTVect3MomXi.Pt() < lCurrentTrck->Pt() ) return;
2076                 // Room for improvement: //FIXME
2077                 // 1. There is a given resolution on pt : maybe release the cut Pt(casc) < Pt(track)*90% ?
2078                 // 2. Apply this cut only when DeltaPhi(casc, track) > 90 deg = when track is in the away-side ?
2079                 // 3. Anti-splitting cut (like in Femto analysis) ?
2080
2081         }// end control loop
2082
2083         // 2nd loop: filling loop
2084         for(Int_t TrckIdx = 0; TrckIdx < (InputEvent())->GetNumberOfTracks() ; TrckIdx++ )
2085         {// Loop over all the tracks of the event
2086
2087                 AliVTrack *lCurrentTrck = dynamic_cast<AliVTrack*>( (InputEvent())->GetTrack( TrckIdx ) );
2088                         if (!lCurrentTrck ) {
2089                                 Printf("ERROR Correl. Study : Could not retrieve a track while looping over the event tracks ...");
2090                                 continue;
2091                         }
2092
2093                 // Room for improvement: //FIXME
2094                 // 1. Loop only on primary tracks ?
2095                 // 2. Exclude the tracks that build the condisdered cascade = the bachelor + the V0 dghters
2096                 //     This may bias the outcome, especially for low multplicity events.
2097                 // Note : For ESD event, track ID == track index.
2098                         if(lCurrentTrck->GetID() == lArrTrackID[0]) continue;
2099                         if(lCurrentTrck->GetID() == lArrTrackID[1]) continue;
2100                         if(lCurrentTrck->GetID() == lArrTrackID[2]) continue;
2101
2102                 TVector3 lTVect3MomTrck(lCurrentTrck->Px(), lCurrentTrck->Py(), lCurrentTrck->Pz() );
2103
2104                 // 2 hypotheses made here :
2105                 //   - The Xi trajectory is a straight line,
2106                 //   - The Xi doesn't loose any energy by crossing the first layer(s) of ITS, if ever;
2107                 //      So, meaning hyp: vect p(Xi) at the emission = vect p(Xi) at the decay vertex
2108                 //      By doing this, we introduce a systematic error on the cascade Phi ...
2109                 // Room for improvement: take into account the curvature of the Xi trajectory //FIXME
2110
2111                 Double_t lHnSpFillVar[5] = {0.};
2112                 lHnSpFillVar[0] = lTVect3MomXi.DeltaPhi(lTVect3MomTrck) * 180.0/TMath::Pi(); // Delta phi(Casc,Track) (deg)
2113                 if(lHnSpFillVar[0] < -50.0) lHnSpFillVar[0] += 360.0;
2114                 lHnSpFillVar[1] = lEtaXi - lCurrentTrck->Eta();                            // Delta eta(Casc,Track)
2115                 lHnSpFillVar[2] = lTVect3MomXi.Pt();                                       // Pt_{Casc}
2116                 lHnSpFillVar[3] = lCurrentTrck->Pt();                                      // Pt_{any track}
2117                 lHnSpFillVar[4] = lInvMassCascade;                                         // Eff. Inv Mass (control var)
2118
2119                 if(      lStrCascType.Contains("Xi-") )      fHnSpAngularCorrXiMinus    ->Fill( lHnSpFillVar );
2120                 else if( lStrCascType.Contains("Xi+") )      fHnSpAngularCorrXiPlus     ->Fill( lHnSpFillVar );
2121                 else if( lStrCascType.Contains("Omega-") )   fHnSpAngularCorrOmegaMinus ->Fill( lHnSpFillVar );
2122                 else if( lStrCascType.Contains("Omega+") )   fHnSpAngularCorrOmegaPlus  ->Fill( lHnSpFillVar );
2123
2124         }// end - Loop over all the tracks in the event
2125
2126 }
2127
2128
2129
2130
2131
2132
2133
2134 //________________________________________________________________________
2135 void AliAnalysisTaskCheckCascade::Terminate(Option_t *)
2136 {
2137   // Draw result to the screen
2138   // Called once at the end of the query
2139
2140   TList *cRetrievedList = 0x0;
2141          cRetrievedList = (TList*)GetOutputData(1);
2142         if(!cRetrievedList){
2143                 Printf("ERROR - AliAnalysisTaskCheckCascade: ouput data container list not available\n"); return;
2144         }
2145
2146   fHistTrackMultiplicity = dynamic_cast<TH1F*> (   cRetrievedList->FindObject("fHistTrackMultiplicity") );
2147   if (!fHistTrackMultiplicity) {
2148                 Printf("ERROR - AliAnalysisTaskCheckCascade: fHistTrackMultiplicity not available\n"); return;
2149         }
2150
2151   fHistCascadeMultiplicity = dynamic_cast<TH1F*> ( cRetrievedList->FindObject("fHistCascadeMultiplicity") );
2152         if (!fHistCascadeMultiplicity) {
2153                 Printf("ERROR - AliAnalysisTaskCheckCascade: fHistCascadeMultiplicity not available\n"); return;
2154         }
2155
2156   fHistMassXiMinus    = dynamic_cast<TH1F*> ( cRetrievedList->FindObject("fHistMassXiMinus") );
2157         if (!fHistMassXiMinus) {
2158                 Printf("ERROR - AliAnalysisTaskCheckCascade: fHistMassXiMinus not available\n"); return;
2159         }
2160   fHistMassXiPlus     = dynamic_cast<TH1F*> ( cRetrievedList->FindObject("fHistMassXiPlus") );
2161         if (!fHistMassXiPlus) {
2162                 Printf("ERROR - AliAnalysisTaskCheckCascade: fHistMassXiPlus not available\n"); return;
2163         }
2164   fHistMassOmegaMinus = dynamic_cast<TH1F*> ( cRetrievedList->FindObject("fHistMassOmegaMinus") );
2165         if (!fHistMassOmegaMinus) {
2166                 Printf("ERROR - AliAnalysisTaskCheckCascade: fHistMassOmegaMinus not available\n"); return;
2167         }
2168   fHistMassOmegaPlus  = dynamic_cast<TH1F*> ( cRetrievedList->FindObject("fHistMassOmegaPlus") );
2169         if (!fHistMassOmegaPlus) {
2170                 Printf("ERROR - AliAnalysisTaskCheckCascade: fHistMassOmegaPlus not available\n"); return;
2171         }
2172
2173   TCanvas *canCheckCascade = new TCanvas("AliAnalysisTaskCheckCascade","CheckCascade overview",10,10,510,510);
2174   canCheckCascade->Divide(2,2);
2175
2176   canCheckCascade->cd(1);
2177   canCheckCascade->cd(1)->SetLogy();
2178   fHistTrackMultiplicity->DrawCopy("H");
2179
2180   canCheckCascade->cd(2);
2181   canCheckCascade->cd(2)->SetLogy();
2182   fHistCascadeMultiplicity->SetMarkerStyle(kOpenSquare);
2183   fHistCascadeMultiplicity->DrawCopy("E");
2184
2185   canCheckCascade->cd(3);
2186   fHistMassXiMinus->SetMarkerStyle(kFullCircle);
2187   fHistMassXiMinus ->SetMarkerSize(0.5);
2188   fHistMassXiMinus->DrawCopy("E");
2189   fHistMassXiPlus ->SetMarkerStyle(kOpenCircle);
2190   fHistMassXiPlus ->SetMarkerSize(0.5);
2191   fHistMassXiPlus ->DrawCopy("ESAME");
2192
2193   canCheckCascade->cd(4);
2194   fHistMassOmegaMinus->SetMarkerStyle(kFullCircle);
2195   fHistMassOmegaMinus->SetMarkerSize(0.5);
2196   fHistMassOmegaMinus->DrawCopy("E");
2197   fHistMassOmegaPlus ->SetMarkerStyle(kOpenCircle);
2198   fHistMassOmegaPlus ->SetMarkerSize(0.5);
2199   fHistMassOmegaPlus ->DrawCopy("ESAME");
2200
2201 }