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