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 three 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
  22 //              3. Rough azimuthal correlation study (Eta, Phi)
  23 //            Adapted to Cascade : A.Maire Mar2008, antonin.maire@ires.in2p3.fr
  24 //            Modified :           A.Maire Aug2009, antonin.maire@ires.in2p3.fr
  25 //-----------------------------------------------------------------
  26
  27
  28
  29 class TTree;
  30 class TParticle;
  31 class TVector3;
  32
  33 //class AliMCEventHandler;
  34 //class AliMCEvent;
  35 //class AliStack;
  36
  37 class AliESDVertex;
  38 class AliAODVertex;
  39 class AliESDv0;
  40 class AliAODv0;
  41
  42 #include <Riostream.h>
  43
  44 #include "TList.h"
  45 #include "TH1.h"
  46 #include "TH2.h"
  47 #include "TH3.h"
  48 #include "THnSparse.h"
  49 #include "TVector3.h"
  50 #include "TCanvas.h"
  51 #include "TMath.h"
  52
  53
  54 #include "AliLog.h"
  55
  56 #include "AliESDEvent.h"
  57 #include "AliAODEvent.h"
  58 //#include "AliCascadeVertexer.h"
  59
  60 #include "AliESDcascade.h"
  61 #include "AliAODcascade.h"
  62
  63 #include "AliAnalysisTaskCheckCascade.h"
  64
  65 ClassImp(AliAnalysisTaskCheckCascade)
  66
  67
  68
  69 //________________________________________________________________________
  70 AliAnalysisTaskCheckCascade::AliAnalysisTaskCheckCascade()
  71   : AliAnalysisTaskSE(), fAnalysisType("ESD"), fCollidingSystems(0),
  72
  73         // - Cascade part initialisation
  74     fListHistCascade(0),
  75     fHistTrackMultiplicity(0), fHistCascadeMultiplicity(0),
  76     fHistVtxStatus(0),
  77
  78     fHistPosTrkgPrimaryVtxX(0), fHistPosTrkgPrimaryVtxY(0), fHistPosTrkgPrimaryVtxZ(0), fHistTrkgPrimaryVtxRadius(0),
  79     fHistPosBestPrimaryVtxX(0), fHistPosBestPrimaryVtxY(0), fHistPosBestPrimaryVtxZ(0), fHistBestPrimaryVtxRadius(0),
  80     f2dHistTrkgPrimVtxVsBestPrimVtx(0),
  81
  82     fHistEffMassXi(0),  fHistChi2Xi(0),
  83     fHistDcaXiDaughters(0), fHistDcaBachToPrimVertex(0), fHistXiCosineOfPointingAngle(0), fHistXiRadius(0),
  84
  85     fHistMassLambdaAsCascDghter(0),
  86     fHistV0Chi2Xi(0),
  87     fHistDcaV0DaughtersXi(0),
  88     fHistDcaV0ToPrimVertexXi(0),
  89     fHistV0CosineOfPointingAngleXi(0),
  90     fHistV0RadiusXi(0),
  91     fHistDcaPosToPrimVertexXi(0), fHistDcaNegToPrimVertexXi(0),
  92
  93     fHistMassXiMinus(0), fHistMassXiPlus(0),
  94     fHistMassOmegaMinus(0), fHistMassOmegaPlus(0),
  95     fHistMassWithCombPIDXiMinus(0), fHistMassWithCombPIDXiPlus(0),
  96     fHistMassWithCombPIDOmegaMinus(0), fHistMassWithCombPIDOmegaPlus(0),
  97
  98     fHistXiTransvMom(0),    fHistXiTotMom(0),
  99     fHistBachTransvMom(0),   fHistBachTotMom(0),
 100
 101     fHistChargeXi(0),
 102     fHistV0toXiCosineOfPointingAngle(0),
 103
 104     fHistRapXi(0), fHistRapOmega(0), fHistEta(0),
 105     fHistTheta(0), fHistPhi(0),
 106
 107     f2dHistArmenteros(0),
 108     f2dHistEffMassLambdaVsEffMassXiMinus(0), f2dHistEffMassXiVsEffMassOmegaMinus(0),
 109     f2dHistEffMassLambdaVsEffMassXiPlus(0), f2dHistEffMassXiVsEffMassOmegaPlus(0),
 110     f2dHistXiRadiusVsEffMassXiMinus(0), f2dHistXiRadiusVsEffMassXiPlus(0),
 111     f2dHistXiRadiusVsEffMassOmegaMinus(0), f2dHistXiRadiusVsEffMassOmegaPlus(0),
 112
 113     f3dHistXiPtVsEffMassVsYXiMinus(0), f3dHistXiPtVsEffMassVsYXiPlus(0),
 114     f3dHistXiPtVsEffMassVsYOmegaMinus(0), f3dHistXiPtVsEffMassVsYOmegaPlus(0),
 115
 116     f3dHistXiPtVsEffMassVsYWithCombPIDXiMinus(0), f3dHistXiPtVsEffMassVsYWithCombPIDXiPlus(0),
 117     f3dHistXiPtVsEffMassVsYWithCombPIDOmegaMinus(0),  f3dHistXiPtVsEffMassVsYWithCombPIDOmegaPlus(0),
 118     f3dHistXiPtVsEffMassVsYWith2CombPIDOmegaMinus(0), f3dHistXiPtVsEffMassVsYWith2CombPIDOmegaPlus(0),
 119
 120     fHnSpAngularCorrXiMinus(0), fHnSpAngularCorrXiPlus(0),
 121     fHnSpAngularCorrOmegaMinus(0), fHnSpAngularCorrOmegaPlus(0)
 122
 123 {
 124   // Dummy Constructor
 125 }
 126
 127
 128
 129
 130
 131
 132
 133
 134 //________________________________________________________________________
 135 AliAnalysisTaskCheckCascade::AliAnalysisTaskCheckCascade(const char *name)
 136   : AliAnalysisTaskSE(name), fAnalysisType("ESD"), fCollidingSystems(0),
 137
 138         // - Cascade part initialisation
 139     fListHistCascade(0),
 140     fHistTrackMultiplicity(0), fHistCascadeMultiplicity(0),
 141     fHistVtxStatus(0),
 142
 143     fHistPosTrkgPrimaryVtxX(0), fHistPosTrkgPrimaryVtxY(0), fHistPosTrkgPrimaryVtxZ(0), fHistTrkgPrimaryVtxRadius(0),
 144     fHistPosBestPrimaryVtxX(0), fHistPosBestPrimaryVtxY(0), fHistPosBestPrimaryVtxZ(0), fHistBestPrimaryVtxRadius(0),
 145     f2dHistTrkgPrimVtxVsBestPrimVtx(0),
 146
 147     fHistEffMassXi(0),  fHistChi2Xi(0),
 148     fHistDcaXiDaughters(0), fHistDcaBachToPrimVertex(0), fHistXiCosineOfPointingAngle(0), fHistXiRadius(0),
 149
 150     fHistMassLambdaAsCascDghter(0),
 151     fHistV0Chi2Xi(0),
 152     fHistDcaV0DaughtersXi(0),
 153     fHistDcaV0ToPrimVertexXi(0),
 154     fHistV0CosineOfPointingAngleXi(0),
 155     fHistV0RadiusXi(0),
 156     fHistDcaPosToPrimVertexXi(0), fHistDcaNegToPrimVertexXi(0),
 157
 158     fHistMassXiMinus(0), fHistMassXiPlus(0),
 159     fHistMassOmegaMinus(0), fHistMassOmegaPlus(0),
 160     fHistMassWithCombPIDXiMinus(0), fHistMassWithCombPIDXiPlus(0),
 161     fHistMassWithCombPIDOmegaMinus(0), fHistMassWithCombPIDOmegaPlus(0),
 162
 163     fHistXiTransvMom(0),    fHistXiTotMom(0),
 164     fHistBachTransvMom(0),   fHistBachTotMom(0),
 165
 166     fHistChargeXi(0),
 167     fHistV0toXiCosineOfPointingAngle(0),
 168
 169     fHistRapXi(0), fHistRapOmega(0), fHistEta(0),
 170     fHistTheta(0), fHistPhi(0),
 171
 172     f2dHistArmenteros(0),
 173     f2dHistEffMassLambdaVsEffMassXiMinus(0), f2dHistEffMassXiVsEffMassOmegaMinus(0),
 174     f2dHistEffMassLambdaVsEffMassXiPlus(0), f2dHistEffMassXiVsEffMassOmegaPlus(0),
 175     f2dHistXiRadiusVsEffMassXiMinus(0), f2dHistXiRadiusVsEffMassXiPlus(0),
 176     f2dHistXiRadiusVsEffMassOmegaMinus(0), f2dHistXiRadiusVsEffMassOmegaPlus(0),
 177
 178     f3dHistXiPtVsEffMassVsYXiMinus(0), f3dHistXiPtVsEffMassVsYXiPlus(0),
 179     f3dHistXiPtVsEffMassVsYOmegaMinus(0), f3dHistXiPtVsEffMassVsYOmegaPlus(0),
 180
 181     f3dHistXiPtVsEffMassVsYWithCombPIDXiMinus(0), f3dHistXiPtVsEffMassVsYWithCombPIDXiPlus(0),
 182     f3dHistXiPtVsEffMassVsYWithCombPIDOmegaMinus(0),  f3dHistXiPtVsEffMassVsYWithCombPIDOmegaPlus(0),
 183     f3dHistXiPtVsEffMassVsYWith2CombPIDOmegaMinus(0), f3dHistXiPtVsEffMassVsYWith2CombPIDOmegaPlus(0),
 184
 185
 186     fHnSpAngularCorrXiMinus(0), fHnSpAngularCorrXiPlus(0),
 187     fHnSpAngularCorrOmegaMinus(0), fHnSpAngularCorrOmegaPlus(0)
 188
 189 {
 190   // Constructor
 191
 192   // Define input and output slots here
 193   // Input slot #0 works with a TChain
 194
 195   // Output slot #0 writes into a TList container (Cascade)
 196   DefineOutput(1, TList::Class());
 197 }
 198
 199
 200
 201
 202
 203
 204 //________________________________________________________________________
 205 void AliAnalysisTaskCheckCascade::UserCreateOutputObjects()
 206 {
 207   // Create histograms
 208   // Called once
 209
 210
 211
 212  fListHistCascade = new TList();
 213
 214
 215         // - General histos
 216
 217 if(! fHistTrackMultiplicity) {
 218         if(fCollidingSystems)// AA collisions
 219                 fHistTrackMultiplicity = new TH1F("fHistTrackMultiplicity",
 220                         "Multiplicity distribution;Number of tracks;Events",
 221                         200, 0, 40000);
 222         else // pp collisions
 223                 fHistTrackMultiplicity = new TH1F("fHistTrackMultiplicity",
 224                         "Track Multiplicity;Nbr of tracks/Evt;Events",
 225                         200, 0, 200);
 226         fListHistCascade->Add(fHistTrackMultiplicity);
 227 }
 228
 229 if(! fHistCascadeMultiplicity) {
 230         if(fCollidingSystems)// AA collisions
 231                 fHistCascadeMultiplicity = new TH1F("fHistCascadeMultiplicity",
 232                         "Multiplicity distribution;Number of Cascades;Events",
 233                         25, 0, 25);
 234         else // pp collisions
 235                 fHistCascadeMultiplicity = new TH1F("fHistCascadeMultiplicity",
 236                         "Cascades per event;Nbr of Cascades/Evt;Events",
 237                         10, 0, 10);
 238         fListHistCascade->Add(fHistCascadeMultiplicity);
 239 }
 240
 241
 242
 243 if(! fHistVtxStatus ){
 244         fHistVtxStatus   = new TH1F( "fHistVtxStatus" , "Does a Trckg Prim.vtx exist ?; true=1 or false=0; Nb of Events" , 4, -1.0, 3.0 );
 245         fListHistCascade->Add(fHistVtxStatus);
 246 }
 247
 248
 249
 250
 251
 252
 253         // - Vertex Positions
 254
 255 if(! fHistPosTrkgPrimaryVtxX ){
 256         fHistPosTrkgPrimaryVtxX   = new TH1F( "fHistPosTrkgPrimaryVtxX" , "Trkg Prim. Vertex Position in x; x (cm); Events" , 200, -0.5, 0.5 );
 257         fListHistCascade->Add(fHistPosTrkgPrimaryVtxX);
 258 }
 259
 260
 261 if(! fHistPosTrkgPrimaryVtxY){
 262         fHistPosTrkgPrimaryVtxY   = new TH1F( "fHistPosTrkgPrimaryVtxY" , "Trkg Prim. Vertex Position in y; y (cm); Events" , 200, -0.5, 0.5 );
 263         fListHistCascade->Add(fHistPosTrkgPrimaryVtxY);
 264 }
 265
 266 if(! fHistPosTrkgPrimaryVtxZ ){
 267         fHistPosTrkgPrimaryVtxZ   = new TH1F( "fHistPosTrkgPrimaryVtxZ" , "Trkg Prim. Vertex Position in z; z (cm); Events" , 100, -15.0, 15.0 );
 268         fListHistCascade->Add(fHistPosTrkgPrimaryVtxZ);
 269 }
 270
 271 if(! fHistTrkgPrimaryVtxRadius ){
 272         fHistTrkgPrimaryVtxRadius  = new TH1F( "fHistTrkgPrimaryVtxRadius",  "Trkg Prim. Vertex radius; r (cm); Events" , 150, 0., 15.0 );
 273         fListHistCascade->Add(fHistTrkgPrimaryVtxRadius);
 274 }
 275
 276
 277
 278
 279 if(! fHistPosBestPrimaryVtxX ){
 280         fHistPosBestPrimaryVtxX   = new TH1F( "fHistPosBestPrimaryVtxX" , "Best Prim. Vertex Position in x; x (cm); Events" , 200, -0.5, 0.5 );
 281         fListHistCascade->Add(fHistPosBestPrimaryVtxX);
 282 }
 283
 284 if(! fHistPosBestPrimaryVtxY){
 285         fHistPosBestPrimaryVtxY   = new TH1F( "fHistPosBestPrimaryVtxY" , "Best Prim. Vertex Position in y; y (cm); Events" , 200, -0.5, 0.5 );
 286         fListHistCascade->Add(fHistPosBestPrimaryVtxY);
 287 }
 288
 289 if(! fHistPosBestPrimaryVtxZ ){
 290         fHistPosBestPrimaryVtxZ   = new TH1F( "fHistPosBestPrimaryVtxZ" , "Best Prim. Vertex Position in z; z (cm); Events" , 100, -15.0, 15.0 );
 291         fListHistCascade->Add(fHistPosBestPrimaryVtxZ);
 292 }
 293
 294 if(! fHistBestPrimaryVtxRadius ){
 295         fHistBestPrimaryVtxRadius  = new TH1F( "fHistBestPrimaryVtxRadius",  "Best Prim.  vertex radius; r (cm); Events" , 150, 0., 15.0 );
 296         fListHistCascade->Add(fHistBestPrimaryVtxRadius);
 297 }
 298
 299 if(! f2dHistTrkgPrimVtxVsBestPrimVtx) {
 300         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.);
 301         fListHistCascade->Add(f2dHistTrkgPrimVtxVsBestPrimVtx);
 302 }
 303
 304
 305
 306
 307 // - Typical histos for cascades
 308
 309
 310 if(! fHistEffMassXi) {
 311      fHistEffMassXi = new TH1F("fHistEffMassXi", "Cascade candidates ; Invariant Mass (GeV/c^{2}) ; Counts", 200, 1.2, 2.0);
 312      fListHistCascade->Add(fHistEffMassXi);
 313 }
 314
 315 if(! fHistChi2Xi ){
 316         fHistChi2Xi = new TH1F("fHistChi2Xi", "Cascade #chi^{2}; #chi^{2}; Number of Cascades", 160, 0, 40);
 317         fListHistCascade->Add(fHistChi2Xi);
 318 }
 319
 320 if(! fHistDcaXiDaughters ){
 321         fHistDcaXiDaughters = new TH1F( "fHistDcaXiDaughters",  "DCA between Xi Daughters; DCA (cm) ; Number of Cascades", 100, 0., 0.5);
 322         fListHistCascade->Add(fHistDcaXiDaughters);
 323 }
 324
 325 if(! fHistDcaBachToPrimVertex) {
 326         fHistDcaBachToPrimVertex = new TH1F("fHistDcaBachToPrimVertex", "DCA of Bach. to Prim. Vertex;DCA (cm);Number of Cascades", 250, 0., 0.25);
 327         fListHistCascade->Add(fHistDcaBachToPrimVertex);
 328 }
 329
 330 if(! fHistXiCosineOfPointingAngle) {
 331         fHistXiCosineOfPointingAngle = new TH1F("fHistXiCosineOfPointingAngle", "Cosine of Xi Pointing Angle; Cos (Xi Point.Angl);Number of Xis", 200, 0.98, 1.0);
 332         fListHistCascade->Add(fHistXiCosineOfPointingAngle);
 333 }
 334
 335 if(! fHistXiRadius ){
 336         fHistXiRadius  = new TH1F( "fHistXiRadius",  "Casc. decay transv. radius; r (cm); Counts" , 200, 0., 20.0 );
 337         fListHistCascade->Add(fHistXiRadius);
 338 }
 339
 340
 341 // - Histos about ~ the "V0 part" of the cascade,  coming by inheritance from AliESDv0
 342
 343
 344
 345 if (! fHistMassLambdaAsCascDghter) {
 346      fHistMassLambdaAsCascDghter = new TH1F("fHistMassLambdaAsCascDghter","#Lambda associated to Casc. candidates;Eff. Mass (GeV/c^{2});Counts", 300,1.00,1.3);
 347     fListHistCascade->Add(fHistMassLambdaAsCascDghter);
 348 }
 349
 350 if (! fHistV0Chi2Xi) {
 351         fHistV0Chi2Xi = new TH1F("fHistV0Chi2Xi", "V0 #chi^{2}, in cascade; #chi^{2};Counts", 160, 0, 40);
 352         fListHistCascade->Add(fHistV0Chi2Xi);
 353 }
 354
 355 if (! fHistDcaV0DaughtersXi) {
 356         fHistDcaV0DaughtersXi = new TH1F("fHistDcaV0DaughtersXi", "DCA between V0 daughters, in cascade;DCA (cm);Number of V0s", 120, 0., 0.6);
 357         fListHistCascade->Add(fHistDcaV0DaughtersXi);
 358 }
 359
 360 if (! fHistDcaV0ToPrimVertexXi) {
 361         fHistDcaV0ToPrimVertexXi = new TH1F("fHistDcaV0ToPrimVertexXi", "DCA of V0 to Prim. Vertex, in cascade;DCA (cm);Number of Cascades", 200, 0., 1.);
 362         fListHistCascade->Add(fHistDcaV0ToPrimVertexXi);
 363 }
 364
 365 if (! fHistV0CosineOfPointingAngleXi) {
 366         fHistV0CosineOfPointingAngleXi = new TH1F("fHistV0CosineOfPointingAngleXi", "Cosine of V0 Pointing Angle, in cascade;Cos(V0 Point. Angl); Counts", 200, 0.98, 1.0);
 367         fListHistCascade->Add(fHistV0CosineOfPointingAngleXi);
 368 }
 369
 370 if (! fHistV0RadiusXi) {
 371         fHistV0RadiusXi = new TH1F("fHistV0RadiusXi", "V0 decay radius, in cascade; radius (cm); Counts", 200, 0, 20);
 372         fListHistCascade->Add(fHistV0RadiusXi);
 373 }
 374
 375 if (! fHistDcaPosToPrimVertexXi) {
 376         fHistDcaPosToPrimVertexXi = new TH1F("fHistDcaPosToPrimVertexXi", "DCA of V0 pos daughter to Prim. Vertex;DCA (cm);Counts", 300, 0, 3);
 377         fListHistCascade->Add(fHistDcaPosToPrimVertexXi);
 378 }
 379
 380 if (! fHistDcaNegToPrimVertexXi) {
 381         fHistDcaNegToPrimVertexXi = new TH1F("fHistDcaNegToPrimVertexXi", "DCA of V0 neg daughter to Prim. Vertex;DCA (cm);Counts", 300, 0, 3);
 382         fListHistCascade->Add(fHistDcaNegToPrimVertexXi);
 383 }
 384
 385
 386
 387
 388         // - Effective mass histos for cascades.
 389 // By cascade hyp
 390 if (! fHistMassXiMinus) {
 391     fHistMassXiMinus = new TH1F("fHistMassXiMinus","#Xi^{-} candidates;M( #Lambda , #pi^{-} ) (GeV/c^{2});Counts", 200,1.2,2.0);
 392     fListHistCascade->Add(fHistMassXiMinus);
 393 }
 394
 395 if (! fHistMassXiPlus) {
 396     fHistMassXiPlus = new TH1F("fHistMassXiPlus","#Xi^{+} candidates;M( #bar{#Lambda}^{0} , #pi^{+} ) (GeV/c^{2});Counts",200,1.2,2.0);
 397     fListHistCascade->Add(fHistMassXiPlus);
 398 }
 399
 400 if (! fHistMassOmegaMinus) {
 401         fHistMassOmegaMinus = new TH1F("fHistMassOmegaMinus","#Omega^{-} candidates;M( #Lambda , K^{-} ) (GeV/c^{2});Counts", 250,1.5,2.5);
 402     fListHistCascade->Add(fHistMassOmegaMinus);
 403 }
 404
 405 if (! fHistMassOmegaPlus) {
 406         fHistMassOmegaPlus = new TH1F("fHistMassOmegaPlus","#Omega^{+} candidates;M( #bar{#Lambda}^{0} , K^{+} ) (GeV/c^{2});Counts", 250,1.5,2.5);
 407     fListHistCascade->Add(fHistMassOmegaPlus);
 408 }
 409
 410 // By cascade hyp + bachelor PID
 411 if (! fHistMassWithCombPIDXiMinus) {
 412     fHistMassWithCombPIDXiMinus = new TH1F("fHistMassWithCombPIDXiMinus","#Xi^{-} candidates, with Bach. comb. PID;M( #Lambda , #pi^{-} ) (GeV/c^{2});Counts", 200,1.2,2.0);
 413     fListHistCascade->Add(fHistMassWithCombPIDXiMinus);
 414 }
 415
 416 if (! fHistMassWithCombPIDXiPlus) {
 417     fHistMassWithCombPIDXiPlus = new TH1F("fHistMassWithCombPIDXiPlus","#Xi^{+} candidates, with Bach. comb. PID;M( #bar{#Lambda}^{0} , #pi^{+} ) (GeV/c^{2});Counts",200,1.2,2.0);
 418     fListHistCascade->Add(fHistMassWithCombPIDXiPlus);
 419 }
 420
 421 if (! fHistMassWithCombPIDOmegaMinus) {
 422         fHistMassWithCombPIDOmegaMinus = new TH1F("fHistMassWithCombPIDOmegaMinus","#Omega^{-} candidates, with Bach. comb. PID;M( #Lambda , K^{-} ) (GeV/c^{2});Counts", 250,1.5,2.5);
 423     fListHistCascade->Add(fHistMassWithCombPIDOmegaMinus);
 424 }
 425
 426 if (! fHistMassWithCombPIDOmegaPlus) {
 427         fHistMassWithCombPIDOmegaPlus = new TH1F("fHistMassWithCombPIDOmegaPlus","#Omega^{+} candidates, with Bach. comb. PID;M( #bar{#Lambda}^{0} , K^{+} ) (GeV/c^{2});Counts", 250,1.5,2.5);
 428     fListHistCascade->Add(fHistMassWithCombPIDOmegaPlus);
 429 }
 430
 431
 432
 433         // - Complements for QA
 434
 435 if(! fHistXiTransvMom ){
 436         fHistXiTransvMom  = new TH1F( "fHistXiTransvMom" , "Xi transverse momentum ; p_{t}(#Xi) (GeV/c); Counts", 100, 0.0, 10.0);
 437         fListHistCascade->Add(fHistXiTransvMom);
 438 }
 439
 440 if(! fHistXiTotMom ){
 441         fHistXiTotMom  = new TH1F( "fHistXiTotMom" , "Xi momentum norm; p_{tot}(#Xi) (GeV/c); Counts", 150, 0.0, 15.0);
 442         fListHistCascade->Add(fHistXiTotMom);
 443 }
 444
 445
 446 if(! fHistBachTransvMom ){
 447         fHistBachTransvMom  = new TH1F( "fHistBachTransvMom" , "Bach. transverse momentum ; p_{t}(Bach.) (GeV/c); Counts", 100, 0.0, 5.0);
 448         fListHistCascade->Add(fHistBachTransvMom);
 449 }
 450
 451 if(! fHistBachTotMom ){
 452         fHistBachTotMom  = new TH1F( "fHistBachTotMom" , "Bach. momentum norm; p_{tot}(Bach.) (GeV/c); Counts", 100, 0.0, 5.0);
 453         fListHistCascade->Add(fHistBachTotMom);
 454 }
 455
 456
 457 if(! fHistChargeXi ){
 458         fHistChargeXi  = new TH1F( "fHistChargeXi" , "Charge of casc. candidates ; Sign ; Counts", 5, -2.0, 3.0);
 459         fListHistCascade->Add(fHistChargeXi);
 460 }
 461
 462
 463 if (! fHistV0toXiCosineOfPointingAngle) {
 464         fHistV0toXiCosineOfPointingAngle = new TH1F("fHistV0toXiCosineOfPointingAngle", "Cos. of V0 Ptng Angl / Xi vtx ;Cos(V0 Point. Angl / Xi vtx); Counts", 100, 0.99, 1.0);
 465         fListHistCascade->Add(fHistV0toXiCosineOfPointingAngle);
 466 }
 467
 468
 469 if(! fHistRapXi ){
 470         fHistRapXi  = new TH1F( "fHistRapXi" , "Rapidity of Xi candidates ; y ; Counts", 200, -5.0, 5.0);
 471         fListHistCascade->Add(fHistRapXi);
 472 }
 473
 474 if(! fHistRapOmega ){
 475         fHistRapOmega  = new TH1F( "fHistRapOmega" , "Rapidity of Omega candidates ; y ; Counts", 200, -5.0, 5.0);
 476         fListHistCascade->Add(fHistRapOmega);
 477 }
 478
 479 if(! fHistEta ){
 480         fHistEta  = new TH1F( "fHistEta" , "Pseudo-rap. of casc. candidates ; #eta ; Counts", 120, -3.0, 3.0);
 481         fListHistCascade->Add(fHistEta);
 482 }
 483
 484 if(! fHistTheta ){
 485         fHistTheta  = new TH1F( "fHistTheta" , "#theta of casc. candidates ; #theta (deg) ; Counts", 180, 0., 180.0);
 486         fListHistCascade->Add(fHistTheta);
 487 }
 488
 489 if(! fHistPhi ){
 490         fHistPhi  = new TH1F( "fHistPhi" , "#phi of casc. candidates ; #phi (deg) ; Counts", 360, 0., 360.);
 491         fListHistCascade->Add(fHistPhi);
 492 }
 493
 494
 495 if(! f2dHistArmenteros) {
 496         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);
 497         fListHistCascade->Add(f2dHistArmenteros);
 498 }
 499
 500 //-------
 501
 502 if(! f2dHistEffMassLambdaVsEffMassXiMinus) {
 503         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);
 504         fListHistCascade->Add(f2dHistEffMassLambdaVsEffMassXiMinus);
 505 }
 506
 507 if(! f2dHistEffMassXiVsEffMassOmegaMinus) {
 508         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);
 509         fListHistCascade->Add(f2dHistEffMassXiVsEffMassOmegaMinus);
 510 }
 511
 512 if(! f2dHistEffMassLambdaVsEffMassXiPlus) {
 513         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);
 514         fListHistCascade->Add(f2dHistEffMassLambdaVsEffMassXiPlus);
 515 }
 516
 517 if(! f2dHistEffMassXiVsEffMassOmegaPlus) {
 518         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);
 519         fListHistCascade->Add(f2dHistEffMassXiVsEffMassOmegaPlus);
 520 }
 521
 522 //-------
 523
 524 if(! f2dHistXiRadiusVsEffMassXiMinus) {
 525         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);
 526         fListHistCascade->Add(f2dHistXiRadiusVsEffMassXiMinus);
 527 }
 528
 529 if(! f2dHistXiRadiusVsEffMassXiPlus) {
 530         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);
 531         fListHistCascade->Add(f2dHistXiRadiusVsEffMassXiPlus);
 532 }
 533
 534 if(! f2dHistXiRadiusVsEffMassOmegaMinus) {
 535         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);
 536         fListHistCascade->Add(f2dHistXiRadiusVsEffMassOmegaMinus);
 537 }
 538
 539 if(! f2dHistXiRadiusVsEffMassOmegaPlus) {
 540         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);
 541         fListHistCascade->Add(f2dHistXiRadiusVsEffMassOmegaPlus);
 542 }
 543
 544
 545 // Part 2 : Raw material for yield extraction -------
 546
 547 if(! f3dHistXiPtVsEffMassVsYXiMinus) {
 548         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);
 549         fListHistCascade->Add(f3dHistXiPtVsEffMassVsYXiMinus);
 550 }
 551
 552 if(! f3dHistXiPtVsEffMassVsYXiPlus) {
 553         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);
 554         fListHistCascade->Add(f3dHistXiPtVsEffMassVsYXiPlus);
 555 }
 556
 557 if(! f3dHistXiPtVsEffMassVsYOmegaMinus) {
 558         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);
 559         fListHistCascade->Add(f3dHistXiPtVsEffMassVsYOmegaMinus);
 560 }
 561
 562 if(! f3dHistXiPtVsEffMassVsYOmegaPlus) {
 563         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);
 564         fListHistCascade->Add(f3dHistXiPtVsEffMassVsYOmegaPlus);
 565 }
 566
 567 //--
 568 if(! f3dHistXiPtVsEffMassVsYWithCombPIDXiMinus) {
 569         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);
 570         fListHistCascade->Add(f3dHistXiPtVsEffMassVsYWithCombPIDXiMinus);
 571 }
 572
 573 if(! f3dHistXiPtVsEffMassVsYWithCombPIDXiPlus) {
 574         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);
 575         fListHistCascade->Add(f3dHistXiPtVsEffMassVsYWithCombPIDXiPlus);
 576 }
 577
 578 if(! f3dHistXiPtVsEffMassVsYWithCombPIDOmegaMinus) {
 579         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);
 580         fListHistCascade->Add(f3dHistXiPtVsEffMassVsYWithCombPIDOmegaMinus);
 581 }
 582
 583 if(! f3dHistXiPtVsEffMassVsYWithCombPIDOmegaPlus) {
 584         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);
 585         fListHistCascade->Add(f3dHistXiPtVsEffMassVsYWithCombPIDOmegaPlus);
 586 }
 587
 588 //--
 589 if(! f3dHistXiPtVsEffMassVsYWith2CombPIDOmegaMinus) {
 590         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);
 591         fListHistCascade->Add(f3dHistXiPtVsEffMassVsYWith2CombPIDOmegaMinus);
 592 }
 593
 594 if(! f3dHistXiPtVsEffMassVsYWith2CombPIDOmegaPlus) {
 595         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);
 596         fListHistCascade->Add(f3dHistXiPtVsEffMassVsYWith2CombPIDOmegaPlus);
 597 }
 598
 599
 600
 601
 602 // Part 3 : Angular correlation study -------
 603
 604 if(! fHnSpAngularCorrXiMinus){
 605         // Delta Phi(Casc,any trck) Vs Delta Eta(Casc,any trck) Vs Casc Pt Vs Pt of the tracks
 606         // Delta Phi  = 360 bins de -180., 180.
 607         // Delta Eta  = 120 bins de -3.0, 3.0
 608         // Pt Cascade = 100 bins de 0., 10.0,
 609         // Pt track = 150 bins de 0., 15.0
 610
 611    Int_t    bins[5] = { 360, 120, 100, 150, 40};
 612    Double_t xmin[5] = {-50., -3.,  0.,  0., 1.30};
 613    Double_t xmax[5] = { 310., 3., 10., 15., 1.34};
 614    fHnSpAngularCorrXiMinus = new THnSparseF("fHnSpAngularCorrXiMinus", "Angular Correlation for #Xi^{-}:", 5, bins, xmin, xmax);
 615         fHnSpAngularCorrXiMinus->GetAxis(0)->SetTitle(" #Delta#phi(Casc,Track) (deg)");
 616         fHnSpAngularCorrXiMinus->GetAxis(1)->SetTitle(" #Delta#eta(Casc,Track)");
 617         fHnSpAngularCorrXiMinus->GetAxis(2)->SetTitle(" Pt_{Casc} (GeV/c)");
 618         fHnSpAngularCorrXiMinus->GetAxis(3)->SetTitle(" Pt_{any track} (GeV/c)");
 619         fHnSpAngularCorrXiMinus->GetAxis(4)->SetTitle(" Eff. Inv Mass (GeV/c^{2})");
 620         fHnSpAngularCorrXiMinus->Sumw2();
 621    fListHistCascade->Add(fHnSpAngularCorrXiMinus);
 622 }
 623
 624 if(! fHnSpAngularCorrXiPlus){
 625         // Delta Phi(Casc,any trck) Vs Delta Eta(Casc,any trck) Vs Casc Pt Vs Pt of the tracks
 626         // Delta Phi  = 360 bins de -180., 180.
 627         // Delta Eta  = 120 bins de -3.0, 3.0
 628         // Pt Cascade = 100 bins de 0., 10.0,
 629         // Pt track = 150 bins de 0., 15.0
 630    Int_t    bins[5] = { 360, 120, 100, 150, 40};
 631    Double_t xmin[5] = {-50., -3.,  0.,  0., 1.30};
 632    Double_t xmax[5] = { 310., 3., 10., 15., 1.34};
 633    fHnSpAngularCorrXiPlus = new THnSparseF("fHnSpAngularCorrXiPlus", "Angular Correlation for #Xi^{+}:", 5, bins, xmin, xmax);
 634         fHnSpAngularCorrXiPlus->GetAxis(0)->SetTitle(" #Delta#phi(Casc,Track) (deg)");
 635         fHnSpAngularCorrXiPlus->GetAxis(1)->SetTitle(" #Delta#eta(Casc,Track)");
 636         fHnSpAngularCorrXiPlus->GetAxis(2)->SetTitle(" Pt_{Casc} (GeV/c)");
 637         fHnSpAngularCorrXiPlus->GetAxis(3)->SetTitle(" Pt_{any track} (GeV/c)");
 638         fHnSpAngularCorrXiPlus->GetAxis(4)->SetTitle(" Eff. Inv Mass (GeV/c^{2})");
 639         fHnSpAngularCorrXiPlus->Sumw2();
 640    fListHistCascade->Add(fHnSpAngularCorrXiPlus);
 641 }
 642
 643 if(! fHnSpAngularCorrOmegaMinus){
 644         // Delta Phi(Casc,any trck) Vs Delta Eta(Casc,any trck) Vs Casc Pt Vs Pt of the tracks
 645         // Delta Phi  = 360 bins de -180., 180.
 646         // Delta Eta  = 120 bins de -3.0, 3.0
 647         // Pt Cascade = 100 bins de 0., 10.0,
 648         // Pt track = 150 bins de 0., 15.0
 649
 650    Int_t    bins[5] = { 360, 120, 100, 150, 40};
 651    Double_t xmin[5] = {-50., -3.,  0.,  0., 1.65};
 652    Double_t xmax[5] = { 310., 3., 10., 15., 1.69};
 653    fHnSpAngularCorrOmegaMinus = new THnSparseF("fHnSpAngularCorrOmegaMinus", "Angular Correlation for #Omega^{-}:", 5, bins, xmin, xmax);
 654         fHnSpAngularCorrOmegaMinus->GetAxis(0)->SetTitle(" #Delta#phi(Casc,Track) (deg)");
 655         fHnSpAngularCorrOmegaMinus->GetAxis(1)->SetTitle(" #Delta#eta(Casc,Track)");
 656         fHnSpAngularCorrOmegaMinus->GetAxis(2)->SetTitle(" Pt_{Casc} (GeV/c)");
 657         fHnSpAngularCorrOmegaMinus->GetAxis(3)->SetTitle(" Pt_{any track} (GeV/c)");
 658         fHnSpAngularCorrOmegaMinus->GetAxis(4)->SetTitle(" Eff. Inv Mass (GeV/c^{2})");
 659         fHnSpAngularCorrOmegaMinus->Sumw2();
 660    fListHistCascade->Add(fHnSpAngularCorrOmegaMinus);
 661 }
 662
 663 if(! fHnSpAngularCorrOmegaPlus){
 664         // Delta Phi(Casc,any trck) Vs Delta Eta(Casc,any trck) Vs Casc Pt Vs Pt of the tracks
 665         // Delta Phi  = 360 bins de -180., 180.
 666         // Delta Eta  = 120 bins de -3.0, 3.0
 667         // Pt Cascade = 100 bins de 0., 10.0,
 668         // Pt track = 150 bins de 0., 15.0
 669    Int_t    bins[5] = { 360, 120, 100, 150, 40};
 670    Double_t xmin[5] = {-50., -3.,  0.,  0., 1.65};
 671    Double_t xmax[5] = { 310., 3., 10., 15., 1.69};
 672    fHnSpAngularCorrOmegaPlus = new THnSparseF("fHnSpAngularCorrOmegaPlus", "Angular Correlation for #Omega^{+}:", 5, bins, xmin, xmax);
 673         fHnSpAngularCorrOmegaPlus->GetAxis(0)->SetTitle(" #Delta#phi(Casc,Track) (deg)");
 674         fHnSpAngularCorrOmegaPlus->GetAxis(1)->SetTitle(" #Delta#eta(Casc,Track)");
 675         fHnSpAngularCorrOmegaPlus->GetAxis(2)->SetTitle(" Pt_{Casc} (GeV/c)");
 676         fHnSpAngularCorrOmegaPlus->GetAxis(3)->SetTitle(" Pt_{any track} (GeV/c)");
 677         fHnSpAngularCorrOmegaPlus->GetAxis(4)->SetTitle(" Eff. Inv Mass (GeV/c^{2})");
 678         fHnSpAngularCorrOmegaPlus->Sumw2();
 679    fListHistCascade->Add(fHnSpAngularCorrOmegaPlus);
 680 }
 681
 682
 683 }// end UserCreateOutputObjects
 684
 685
 686
 687
 688
 689
 690 //________________________________________________________________________
 691 void AliAnalysisTaskCheckCascade::UserExec(Option_t *)
 692 {
 693   // Main loop
 694   // Called for each event
 695
 696         AliESDEvent *lESDevent = 0x0;
 697         AliAODEvent *lAODevent = 0x0;
 698         Int_t ncascades = -1;
 699
 700
 701   // Connect to the InputEvent
 702   // After these lines, we should have an ESD/AOD event + the number of cascades in it.
 703
 704   if(fAnalysisType == "ESD"){
 705         lESDevent = dynamic_cast<AliESDEvent*>( InputEvent() );
 706         if (!lESDevent) {
 707                 Printf("ERROR: lESDevent not available \n");
 708                 return;
 709         }
 710         ncascades = lESDevent->GetNumberOfCascades();
 711   }
 712
 713   if(fAnalysisType == "AOD"){
 714           lAODevent = dynamic_cast<AliAODEvent*>( InputEvent() );
 715         if (!lAODevent) {
 716                 Printf("ERROR: lAODevent not available \n");
 717                 return;
 718         }
 719         ncascades = lAODevent->GetNumberOfCascades();
 720         // printf("Number of cascade(s) = %d \n", ncascades);
 721   }
 722
 723   // ---------------------------------------------------------------
 724   // I - Initialisation of the local variables that will be needed
 725
 726
 727         // - 0th part of initialisation : around primary vertex ...
 728   Short_t  lStatusTrackingPrimVtx  = -2;
 729   Double_t lTrkgPrimaryVtxPos[3]   = {-100.0, -100.0, -100.0};
 730   Double_t lTrkgPrimaryVtxRadius3D = -500.0;
 731
 732   Double_t lBestPrimaryVtxPos[3]   = {-100.0, -100.0, -100.0};
 733   Double_t lBestPrimaryVtxRadius3D = -500.0;
 734
 735         // - 1st part of initialisation : variables needed to store AliESDCascade data members
 736   Double_t lEffMassXi      = 0. ;
 737   Double_t lChi2Xi         = 0. ;
 738   Double_t lDcaXiDaughters = 0. ;
 739   Double_t lXiCosineOfPointingAngle = 0. ;
 740   Double_t lPosXi[3] = { -1000.0, -1000.0, -1000.0 };
 741   Double_t lXiRadius = 0. ;
 742
 743
 744         // - 2nd part of initialisation : about V0 part in cascades
 745
 746   Double_t lInvMassLambdaAsCascDghter = 0.;
 747   Double_t lV0Chi2Xi         = 0. ;
 748   Double_t lDcaV0DaughtersXi = 0.;
 749
 750   Double_t lDcaBachToPrimVertexXi = 0., lDcaV0ToPrimVertexXi = 0.;
 751   Double_t lDcaPosToPrimVertexXi  = 0.;
 752   Double_t lDcaNegToPrimVertexXi  = 0.;
 753   Double_t lV0CosineOfPointingAngleXi = 0. ;
 754   Double_t lPosV0Xi[3] = { -1000. , -1000., -1000. }; // Position of VO coming from cascade
 755   Double_t lV0RadiusXi = -1000.0;
 756   Double_t lV0quality  = 0.;
 757
 758
 759         // - 3rd part of initialisation : Effective masses
 760
 761   Double_t lInvMassXiMinus    = 0.;
 762   Double_t lInvMassXiPlus     = 0.;
 763   Double_t lInvMassOmegaMinus = 0.;
 764   Double_t lInvMassOmegaPlus  = 0.;
 765
 766         // - 4th part of initialisation : PID treatment
 767   Bool_t   lIsPosInXiProton    = kFALSE;
 768   Bool_t   lIsPosInOmegaProton = kFALSE;
 769   Bool_t   lIsNegInXiProton    = kFALSE;
 770   Bool_t   lIsNegInOmegaProton = kFALSE;
 771   Bool_t   lIsBachelorKaon     = kFALSE;
 772   Bool_t   lIsBachelorPion     = kFALSE;
 773
 774         // - 5th part of initialisation : extra info for QA
 775
 776   Double_t lXiMomX       = 0., lXiMomY = 0., lXiMomZ = 0.;
 777   Double_t lXiTransvMom  = 0. ;
 778   Double_t lXiTotMom     = 0. ;
 779
 780   Double_t lBachMomX       = 0., lBachMomY  = 0., lBachMomZ   = 0.;
 781   Double_t lBachTransvMom  = 0.;
 782   Double_t lBachTotMom     = 0.;
 783
 784   Short_t  lChargeXi = -2;
 785   Double_t lV0toXiCosineOfPointingAngle = 0. ;
 786
 787   Double_t lRapXi   = -20.0, lRapOmega = -20.0,  lEta = -20.0, lTheta = 360., lPhi = 720. ;
 788   Double_t lAlphaXi = -200., lPtArmXi  = -200.0;
 789
 790   TVector3 lTVect3MomXi(0.,0.,0.);
 791   Int_t    lArrTrackID[3] = {-1, -1, -1};
 792
 793
 794
 795   //-------------------------------------------------
 796   // O - Cascade vertexer (ESD)
 797
 798         //   if(fAnalysisType == "ESD" ){
 799         //      lESDevent->ResetCascades();
 800         //      AliCascadeVertexer CascVtxer;
 801         //      CascVtxer.V0sTracks2CascadeVertices(lESDevent);
 802         //     }
 803
 804
 805   // ---------------------------------------------------------------
 806   // I - General histos (filled for any event)
 807
 808   fHistTrackMultiplicity  ->Fill( (InputEvent())->GetNumberOfTracks() );
 809   fHistCascadeMultiplicity->Fill( ncascades );
 810
 811
 812   for (Int_t iXi = 0; iXi < ncascades; iXi++)
 813   {// This is the begining of the Cascade loop (ESD or AOD)
 814
 815
 816   if(fAnalysisType == "ESD"){
 817
 818   // -------------------------------------
 819   // II - Calcultaion Part dedicated to Xi vertices (ESD)
 820
 821         AliESDcascade *xi = lESDevent->GetCascade(iXi);
 822         if (!xi) continue;
 823
 824         // Just to know which file is currently open : locate the file containing Xi
 825         // cout << "Name of the file containing Xi candidate(s) :"
 826         //      << fInputHandler->GetTree()->GetCurrentFile()->GetName()
 827         //      << endl;
 828
 829
 830                 // - II.Step 1 : Characteristics of the event : prim. Vtx + magnetic field (ESD)
 831                 //-------------
 832
 833         // For new code (v4-16-Release-Rev06 or trunk)
 834
 835         const AliESDVertex *lPrimaryTrackingVtx = lESDevent->GetPrimaryVertexTracks();
 836                 // get the vtx stored in ESD found with tracks
 837
 838                 lPrimaryTrackingVtx->GetXYZ( lTrkgPrimaryVtxPos );
 839                 lTrkgPrimaryVtxRadius3D = TMath::Sqrt(  lTrkgPrimaryVtxPos[0] * lTrkgPrimaryVtxPos[0] +
 840                                                         lTrkgPrimaryVtxPos[1] * lTrkgPrimaryVtxPos[1] +
 841                                                         lTrkgPrimaryVtxPos[2] * lTrkgPrimaryVtxPos[2] );
 842
 843                 lStatusTrackingPrimVtx = lPrimaryTrackingVtx->GetStatus();
 844
 845
 846         const AliESDVertex *lPrimaryBestVtx = lESDevent->GetPrimaryVertex();
 847                 // get the best primary vertex available for the event
 848                 // As done in AliCascadeVertexer, we keep the one which is the best one available.
 849                 // between : Tracking vertex > SPD vertex > TPC vertex > default SPD vertex
 850                 // This one will be used for next calculations (DCA essentially)
 851
 852                 lPrimaryBestVtx->GetXYZ( lBestPrimaryVtxPos );
 853                 lBestPrimaryVtxRadius3D = TMath::Sqrt(  lBestPrimaryVtxPos[0] * lBestPrimaryVtxPos[0] +
 854                                                         lBestPrimaryVtxPos[1] * lBestPrimaryVtxPos[1] +
 855                                                         lBestPrimaryVtxPos[2] * lBestPrimaryVtxPos[2] );
 856
 857
 858         // For older evts
 859
 860         //      const AliESDVertex *lPrimaryTrackingVtx = lESDevent->GetPrimaryVertexTracks();
 861         //              get the vtx stored in ESD found with tracks
 862         //      Double_t lTrkgPrimaryVtxPos[3] = {-100.0, -100.0, -100.0};
 863         //              lPrimaryTrackingVtx->GetXYZ( lTrkgPrimaryVtxPos );
 864         //              Double_t lTrkgPrimaryVtxRadius3D = TMath::Sqrt( lTrkgPrimaryVtxPos[0]*lTrkgPrimaryVtxPos[0] +
 865         //                                              lTrkgPrimaryVtxPos[1] * lTrkgPrimaryVtxPos[1] +
 866         //                                              lTrkgPrimaryVtxPos[2] * lTrkgPrimaryVtxPos[2] );
 867         //
 868         //      const AliESDVertex *lPrimarySPDVtx = lESDevent->GetPrimaryVertexSPD();
 869         //              get the vtx found by exclusive use of SPD
 870         //      Double_t lSPDPrimaryVtxPos[3] = {-100.0, -100.0, -100.0};
 871         //              lPrimarySPDVtx->GetXYZ( lSPDPrimaryVtxPos );
 872         //
 873         //      // As done in AliCascadeVertexer, we keep, between both retrieved vertices,
 874         //      // the one which is the best one available.
 875         //      // This one will be used for next calculations (DCA essentially)
 876         //
 877         //              Double_t lBestPrimaryVtxPos[3] = {-100.0, -100.0, -100.0};
 878         //              if( lPrimaryTrackingVtx->GetStatus() ) { // if tracking vtx = ok
 879         //                      lBestPrimaryVtxPos[0] = lTrkgPrimaryVtxPos[0];
 880         //                      lBestPrimaryVtxPos[1] = lTrkgPrimaryVtxPos[1];
 881         //                      lBestPrimaryVtxPos[2] = lTrkgPrimaryVtxPos[2];
 882         //              }
 883         //              else{
 884         //                      lBestPrimaryVtxPos[0] = lSPDPrimaryVtxPos[0];
 885         //                      lBestPrimaryVtxPos[1] = lSPDPrimaryVtxPos[1];
 886         //                      lBestPrimaryVtxPos[2] = lSPDPrimaryVtxPos[2];
 887         //              }
 888         //
 889         //      Double_t lBestPrimaryVtxRadius3D = TMath::Sqrt( lBestPrimaryVtxPos[0]*lBestPrimaryVtxPos[0] +
 890         //                                              lBestPrimaryVtxPos[1] * lBestPrimaryVtxPos[1] +
 891         //                                              lBestPrimaryVtxPos[2] * lBestPrimaryVtxPos[2] );
 892
 893         // Back to normal
 894
 895         Double_t lMagneticField = lESDevent->GetMagneticField( );
 896
 897
 898
 899                 // - II.Step 2 : Assigning the necessary variables for specific AliESDcascade data members (ESD)
 900                 //-------------
 901         lV0quality = 0.;
 902         xi->ChangeMassHypothesis(lV0quality , 3312); // default working hypothesis : cascade = Xi- decay
 903
 904         lEffMassXi                      = xi->GetEffMassXi();
 905         lChi2Xi                         = xi->GetChi2Xi();
 906         lDcaXiDaughters                 = xi->GetDcaXiDaughters();
 907         lXiCosineOfPointingAngle        = xi->GetCascadeCosineOfPointingAngle( lBestPrimaryVtxPos[0],
 908                                                                                 lBestPrimaryVtxPos[1],
 909                                                                                 lBestPrimaryVtxPos[2] );
 910                 // Take care : the best available vertex should be used (like in AliCascadeVertexer)
 911
 912         xi->GetXYZcascade( lPosXi[0],  lPosXi[1], lPosXi[2] );
 913         lXiRadius                       = TMath::Sqrt( lPosXi[0]*lPosXi[0]  +  lPosXi[1]*lPosXi[1] );
 914
 915
 916
 917                 // - II.Step 3 : around the tracks : Bach + V0 (ESD)
 918                 // ~ Necessary variables for ESDcascade data members coming from the ESDv0 part (inheritance)
 919                 //-------------
 920
 921                 UInt_t lIdxPosXi        = (UInt_t) TMath::Abs( xi->GetPindex() );
 922                 UInt_t lIdxNegXi        = (UInt_t) TMath::Abs( xi->GetNindex() );
 923                 UInt_t lBachIdx         = (UInt_t) TMath::Abs( xi->GetBindex() );
 924                         // Care track label can be negative in MC production (linked with the track quality)
 925                         // However = normally, not the case for track index ...
 926
 927         AliESDtrack *pTrackXi           = lESDevent->GetTrack( lIdxPosXi );
 928         AliESDtrack *nTrackXi           = lESDevent->GetTrack( lIdxNegXi );
 929         AliESDtrack *bachTrackXi        = lESDevent->GetTrack( lBachIdx );
 930         if (!pTrackXi || !nTrackXi || !bachTrackXi ) {
 931                 Printf("ERROR: Could not retrieve one of the 3 ESD daughter tracks of the cascade ...");
 932                 continue;
 933         }
 934
 935         lInvMassLambdaAsCascDghter      = xi->GetEffMass();
 936                 // This value shouldn't change, whatever the working hyp. is : Xi-, Xi+, Omega-, Omega+
 937         lDcaV0DaughtersXi               = xi->GetDcaV0Daughters();
 938         lV0Chi2Xi                       = xi->GetChi2V0();
 939
 940         lV0CosineOfPointingAngleXi      = xi->GetV0CosineOfPointingAngle( lBestPrimaryVtxPos[0],
 941                                                                           lBestPrimaryVtxPos[1],
 942                                                                           lBestPrimaryVtxPos[2] );
 943
 944         lDcaV0ToPrimVertexXi            = xi->GetD( lBestPrimaryVtxPos[0],
 945                                                     lBestPrimaryVtxPos[1],
 946                                                     lBestPrimaryVtxPos[2] );
 947
 948         lDcaBachToPrimVertexXi = TMath::Abs( bachTrackXi->GetD( lBestPrimaryVtxPos[0],
 949                                                                 lBestPrimaryVtxPos[1],
 950                                                                 lMagneticField  ) );
 951                                         // Note : AliExternalTrackParam::GetD returns an algebraic value ...
 952
 953                 xi->GetXYZ( lPosV0Xi[0],  lPosV0Xi[1], lPosV0Xi[2] );
 954         lV0RadiusXi             = TMath::Sqrt( lPosV0Xi[0]*lPosV0Xi[0]  +  lPosV0Xi[1]*lPosV0Xi[1] );
 955
 956         lDcaPosToPrimVertexXi   = TMath::Abs( pTrackXi  ->GetD( lBestPrimaryVtxPos[0],
 957                                                                 lBestPrimaryVtxPos[1],
 958                                                                 lMagneticField  )     );
 959
 960         lDcaNegToPrimVertexXi   = TMath::Abs( nTrackXi  ->GetD( lBestPrimaryVtxPos[0],
 961                                                                 lBestPrimaryVtxPos[1],
 962                                                                 lMagneticField  )     );
 963
 964
 965                 // - II.Step 4 : around effective masses (ESD)
 966                 // ~ change mass hypotheses to cover all the possibilities :  Xi-/+, Omega -/+
 967                 //-------------
 968
 969
 970         if( bachTrackXi->Charge() < 0 ) {
 971                 lV0quality = 0.;
 972                 xi->ChangeMassHypothesis(lV0quality , 3312);
 973                         // Calculate the effective mass of the Xi- candidate.
 974                         // pdg code 3312 = Xi-
 975                 lInvMassXiMinus = xi->GetEffMassXi();
 976
 977                 lV0quality = 0.;
 978                 xi->ChangeMassHypothesis(lV0quality , 3334);
 979                         // Calculate the effective mass of the Xi- candidate.
 980                         // pdg code 3334 = Omega-
 981                 lInvMassOmegaMinus = xi->GetEffMassXi();
 982
 983                 lV0quality = 0.;
 984                 xi->ChangeMassHypothesis(lV0quality , 3312);    // Back to default hyp.
 985         }// end if negative bachelor
 986
 987
 988         if( bachTrackXi->Charge() >  0 ){
 989                 lV0quality = 0.;
 990                 xi->ChangeMassHypothesis(lV0quality , -3312);
 991                         // Calculate the effective mass of the Xi+ candidate.
 992                         // pdg code -3312 = Xi+
 993                 lInvMassXiPlus = xi->GetEffMassXi();
 994
 995                 lV0quality = 0.;
 996                 xi->ChangeMassHypothesis(lV0quality , -3334);
 997                         // Calculate the effective mass of the Xi+ candidate.
 998                         // pdg code -3334  = Omega+
 999                 lInvMassOmegaPlus = xi->GetEffMassXi();
1000
1001                 lV0quality = 0.;
1002                 xi->ChangeMassHypothesis(lV0quality , -3312);   // Back to "default" hyp.
1003         }// end if positive bachelor
1004
1005
1006
1007                 // - II.Step 5 : PID on the bachelor
1008                 //-------------
1009
1010         // Reasonable guess for the priors for the cascade track sample
1011         Double_t lPriorsGuessXi[5]    = {0.0, 0.0, 2, 0, 1};
1012         Double_t lPriorsGuessOmega[5] = {0.0, 0.0, 1, 1, 1};
1013
1014         // VO positive daughter PID
1015         AliPID pPidXi;          pPidXi.SetPriors(    lPriorsGuessXi    );
1016         AliPID pPidOmega;       pPidOmega.SetPriors( lPriorsGuessOmega );
1017
1018         if( pTrackXi->IsOn(AliESDtrack::kESDpid) ){  // Combined PID exists
1019                 Double_t r[10] = {0.}; pTrackXi->GetESDpid(r);
1020                 pPidXi.SetProbabilities(r);
1021                 pPidOmega.SetProbabilities(r);
1022                 // Check if the V0 positive track is a proton (case for Xi-)
1023                 Double_t pproton = pPidXi.GetProbability(AliPID::kProton);
1024                 if (pproton > pPidXi.GetProbability(AliPID::kElectron) &&
1025                     pproton > pPidXi.GetProbability(AliPID::kMuon)     &&
1026                     pproton > pPidXi.GetProbability(AliPID::kPion)     &&
1027                     pproton > pPidXi.GetProbability(AliPID::kKaon)     )     lIsPosInXiProton = kTRUE;
1028                 // Check if the V0 positive track is a proton (case for Omega-)
1029                 pproton = 0.;
1030                          pproton = pPidOmega.GetProbability(AliPID::kProton);
1031                 if (pproton > pPidOmega.GetProbability(AliPID::kElectron) &&
1032                     pproton > pPidOmega.GetProbability(AliPID::kMuon)     &&
1033                     pproton > pPidOmega.GetProbability(AliPID::kPion)     &&
1034                     pproton > pPidOmega.GetProbability(AliPID::kKaon)     )  lIsPosInOmegaProton = kTRUE;
1035         }// end if V0 positive track with existing combined PID
1036
1037         // VO negative daughter PID
1038         AliPID nPidXi;          nPidXi.SetPriors(    lPriorsGuessXi    );
1039         AliPID nPidOmega;       nPidOmega.SetPriors( lPriorsGuessOmega );
1040
1041         if( nTrackXi->IsOn(AliESDtrack::kESDpid) ){  // Combined PID exists
1042                 Double_t r[10] = {0.}; nTrackXi->GetESDpid(r);
1043                 nPidXi.SetProbabilities(r);
1044                 nPidOmega.SetProbabilities(r);
1045                 // Check if the V0 negative track is an anti-proton (case for Xi+)
1046                 Double_t pproton = nPidXi.GetProbability(AliPID::kProton);
1047                 if (pproton > nPidXi.GetProbability(AliPID::kElectron) &&
1048                     pproton > nPidXi.GetProbability(AliPID::kMuon)     &&
1049                     pproton > nPidXi.GetProbability(AliPID::kPion)     &&
1050                     pproton > nPidXi.GetProbability(AliPID::kKaon)     )     lIsNegInXiProton = kTRUE;
1051                 // Check if the V0 negative track is an anti-proton (case for Omega+)
1052                 pproton = 0.;
1053                          pproton = nPidOmega.GetProbability(AliPID::kProton);
1054                 if (pproton > nPidOmega.GetProbability(AliPID::kElectron) &&
1055                     pproton > nPidOmega.GetProbability(AliPID::kMuon)     &&
1056                     pproton > nPidOmega.GetProbability(AliPID::kPion)     &&
1057                     pproton > nPidOmega.GetProbability(AliPID::kKaon)     )  lIsNegInOmegaProton = kTRUE;
1058         }// end if V0 negative track with existing combined PID
1059
1060
1061         // Bachelor PID
1062         AliPID bachPidXi;       bachPidXi.SetPriors(    lPriorsGuessXi    );
1063         AliPID bachPidOmega;    bachPidOmega.SetPriors( lPriorsGuessOmega );
1064
1065         if( bachTrackXi->IsOn(AliESDtrack::kESDpid) ){  // Combined PID exists
1066                 Double_t r[10] = {0.}; bachTrackXi->GetESDpid(r);
1067                 bachPidXi.SetProbabilities(r);
1068                 bachPidOmega.SetProbabilities(r);
1069                 // Check if the bachelor track is a pion
1070                 Double_t ppion = bachPidXi.GetProbability(AliPID::kPion);
1071                 if (ppion > bachPidXi.GetProbability(AliPID::kElectron) &&
1072                     ppion > bachPidXi.GetProbability(AliPID::kMuon)     &&
1073                     ppion > bachPidXi.GetProbability(AliPID::kKaon)     &&
1074                     ppion > bachPidXi.GetProbability(AliPID::kProton)   )     lIsBachelorPion = kTRUE;
1075                 // Check if the bachelor track is a kaon
1076                 Double_t pkaon = bachPidOmega.GetProbability(AliPID::kKaon);
1077                 if (pkaon > bachPidOmega.GetProbability(AliPID::kElectron) &&
1078                     pkaon > bachPidOmega.GetProbability(AliPID::kMuon)     &&
1079                     pkaon > bachPidOmega.GetProbability(AliPID::kPion)     &&
1080                     pkaon > bachPidOmega.GetProbability(AliPID::kProton)   )  lIsBachelorKaon = kTRUE;
1081         }// end if bachelor track with existing combined PID
1082
1083
1084
1085                 // - II.Step 6 : extra info for QA (ESD)
1086                 // miscellaneous pieces of info that may help regarding data quality assessment.
1087                 //-------------
1088
1089         xi->GetPxPyPz( lXiMomX, lXiMomY, lXiMomZ );
1090                 lXiTransvMom    = TMath::Sqrt( lXiMomX*lXiMomX   + lXiMomY*lXiMomY );
1091                 lXiTotMom       = TMath::Sqrt( lXiMomX*lXiMomX   + lXiMomY*lXiMomY   + lXiMomZ*lXiMomZ );
1092
1093         xi->GetBPxPyPz(  lBachMomX,  lBachMomY,  lBachMomZ );
1094                 lBachTransvMom  = TMath::Sqrt( lBachMomX*lBachMomX   + lBachMomY*lBachMomY );
1095                 lBachTotMom     = TMath::Sqrt( lBachMomX*lBachMomX   + lBachMomY*lBachMomY  +  lBachMomZ*lBachMomZ  );
1096
1097         lChargeXi = xi->Charge();
1098
1099         lV0toXiCosineOfPointingAngle = xi->GetV0CosineOfPointingAngle( lPosXi[0], lPosXi[1], lPosXi[2] );
1100
1101         lRapXi    = xi->RapXi();
1102         lRapOmega = xi->RapOmega();
1103         lEta      = xi->Eta();
1104         lTheta    = xi->Theta() *180.0/TMath::Pi();
1105         lPhi      = xi->Phi()   *180.0/TMath::Pi();
1106         lAlphaXi  = xi->AlphaXi();
1107         lPtArmXi  = xi->PtArmXi();
1108
1109
1110                 // II.Step 7 - Azimuthal correlation study
1111                 //-------------
1112
1113         lTVect3MomXi.SetXYZ( lXiMomX, lXiMomY, lXiMomZ );
1114         lArrTrackID[0] = pTrackXi   ->GetID();
1115         lArrTrackID[1] = nTrackXi   ->GetID();
1116         lArrTrackID[2] = bachTrackXi->GetID();
1117
1118
1119   }// end of ESD treatment
1120
1121
1122   if(fAnalysisType == "AOD"){
1123
1124         // -------------------------------------
1125         // II - Calcultaion Part dedicated to Xi vertices (ESD)
1126
1127         const AliAODcascade *xi = lAODevent->GetCascade(iXi);
1128         if (!xi) continue;
1129
1130         // Just to know which file is currently open : locate the file containing Xi
1131         // cout << "Name of the file containing Xi candidate(s) :" <<  fesdH->GetTree()->GetCurrentFile()->GetName() << endl;
1132
1133
1134                 // - II.Step 1 : Characteristics of the event : prim. Vtx + magnetic field (AOD)
1135                 //-------------
1136
1137         lTrkgPrimaryVtxPos[0]   = -100.0;
1138         lTrkgPrimaryVtxPos[1]   = -100.0;
1139         lTrkgPrimaryVtxPos[2]   = -100.0;
1140         lTrkgPrimaryVtxRadius3D = -500. ;
1141         // We don't have the different prim. vertex at the AOD level -> nothing to do.
1142
1143         const AliAODVertex *lPrimaryBestVtx = lAODevent->GetPrimaryVertex();
1144         // get the best primary vertex available for the event
1145         // We may keep the one which is the best one available = GetVertex(0)
1146         // Pb with pile-up to expect
1147         // This one will be used for next calculations (DCA essentially)
1148
1149         lPrimaryBestVtx->GetXYZ( lBestPrimaryVtxPos );
1150         lBestPrimaryVtxRadius3D = TMath::Sqrt(  lBestPrimaryVtxPos[0] * lBestPrimaryVtxPos[0] +
1151                                                 lBestPrimaryVtxPos[1] * lBestPrimaryVtxPos[1] +
1152                                                 lBestPrimaryVtxPos[2] * lBestPrimaryVtxPos[2] );
1153
1154
1155                 // - II.Step 2 : Assigning the necessary variables for specific AliAODcascade data members (AOD)
1156                 //-------------
1157
1158         lEffMassXi                      = xi->MassXi(); // default working hypothesis : cascade = Xi- decay
1159         lChi2Xi                         = xi->Chi2Xi();
1160         lDcaXiDaughters                 = xi->DcaXiDaughters();
1161         lXiCosineOfPointingAngle        = xi->CosPointingAngleXi( lBestPrimaryVtxPos[0],
1162                                                                   lBestPrimaryVtxPos[1],
1163                                                                   lBestPrimaryVtxPos[2] );
1164                                         // Take care :
1165                                         // the best available vertex should be used (like in AliCascadeVertexer)
1166
1167                 lPosXi[0] = xi->DecayVertexXiX();
1168                 lPosXi[1] = xi->DecayVertexXiY();
1169                 lPosXi[2] = xi->DecayVertexXiZ();
1170         lXiRadius = TMath::Sqrt( lPosXi[0]*lPosXi[0]  +  lPosXi[1]*lPosXi[1] );
1171
1172
1173                 // - II.Step 3 : around the tracks : Bach + V0 (AOD)
1174                 // ~ Necessary variables for AODcascade data members coming from the AODv0 part (inheritance)
1175                 //-------------
1176
1177         lChargeXi                       = xi->ChargeXi();
1178
1179         if( lChargeXi < 0)
1180           lInvMassLambdaAsCascDghter    = xi->MassLambda();
1181         else
1182           lInvMassLambdaAsCascDghter    = xi->MassAntiLambda();
1183
1184         lDcaV0DaughtersXi               = xi->DcaV0Daughters();
1185         lV0Chi2Xi                       = xi->Chi2V0();
1186
1187         lV0CosineOfPointingAngleXi      = xi->CosPointingAngle( lBestPrimaryVtxPos );
1188         lDcaV0ToPrimVertexXi            = xi->DcaV0ToPrimVertex();
1189
1190         lDcaBachToPrimVertexXi          = xi->DcaBachToPrimVertex();
1191
1192
1193                 lPosV0Xi[0] = xi->DecayVertexV0X();
1194                 lPosV0Xi[1] = xi->DecayVertexV0Y();
1195                 lPosV0Xi[2] = xi->DecayVertexV0Z();
1196         lV0RadiusXi     = TMath::Sqrt( lPosV0Xi[0]*lPosV0Xi[0]  +  lPosV0Xi[1]*lPosV0Xi[1] );
1197
1198         lDcaPosToPrimVertexXi           = xi->DcaPosToPrimVertex();
1199         lDcaNegToPrimVertexXi           = xi->DcaNegToPrimVertex();
1200
1201
1202                 // - II.Step 4 : around effective masses (AOD)
1203                 // ~ change mass hypotheses to cover all the possibilities :  Xi-/+, Omega -/+
1204                 //-------------
1205
1206         if( lChargeXi < 0 )             lInvMassXiMinus         = xi->MassXi();
1207         if( lChargeXi > 0 )             lInvMassXiPlus          = xi->MassXi();
1208         if( lChargeXi < 0 )             lInvMassOmegaMinus      = xi->MassOmega();
1209         if( lChargeXi > 0 )             lInvMassOmegaPlus       = xi->MassOmega();
1210
1211
1212                 // - II.Step 5 : PID on the bachelor (To be developed ...)
1213                 //-------------
1214
1215         /*
1216         // Reasonable guess for the priors for the cascade track sample
1217         Double_t lPriorsGuessXi[5]    = {0.0, 0.0, 2, 0, 1};
1218         Double_t lPriorsGuessOmega[5] = {0.0, 0.0, 1, 1, 1};
1219         AliPID bachPidXi;       bachPidXi.SetPriors(    lPriorsGuessXi    );
1220         AliPID bachPidOmega;    bachPidOmega.SetPriors( lPriorsGuessOmega );
1221
1222         const AliAODTrack *bachTrackXi = lAODevent->GetTrack( xi->GetBachID() ); // FIXME : GetBachID not implemented ?
1223
1224         if( bachTrackXi->IsOn(AliESDtrack::kESDpid) ){  // Combined PID exists, the AOD flags = a copy of the ESD ones
1225                 Double_t r[10]; bachTrackXi->GetPID(r);
1226                 bachPidXi.SetProbabilities(r);
1227                 bachPidOmega.SetProbabilities(r);
1228                 // Check if the bachelor track is a pion
1229                 Double_t ppion = bachPidXi.GetProbability(AliPID::kPion);
1230                 if (ppion > bachPidXi.GetProbability(AliPID::kElectron) &&
1231                     ppion > bachPidXi.GetProbability(AliPID::kMuon)     &&
1232                     ppion > bachPidXi.GetProbability(AliPID::kKaon)     &&
1233                     ppion > bachPidXi.GetProbability(AliPID::kProton)   )     lIsBachelorPion = kTRUE;
1234                 // Check if the bachelor track is a kaon
1235                 Double_t pkaon = bachPidOmega.GetProbability(AliPID::kKaon);
1236                 if (pkaon > bachPidOmega.GetProbability(AliPID::kElectron) &&
1237                     pkaon > bachPidOmega.GetProbability(AliPID::kMuon)     &&
1238                     pkaon > bachPidOmega.GetProbability(AliPID::kPion)     &&
1239                     pkaon > bachPidOmega.GetProbability(AliPID::kProton)   )  lIsBachelorKaon = kTRUE;
1240
1241         }// end if bachelor track with existing combined PID
1242         */
1243
1244
1245                 // - II.Step 6 : extra info for QA (AOD)
1246                 // miscellaneous pieces onf info that may help regarding data quality assessment.
1247                 //-------------
1248
1249                 lXiMomX = xi->MomXiX();
1250                 lXiMomY = xi->MomXiY();
1251                 lXiMomZ = xi->MomXiZ();
1252         lXiTransvMom    = TMath::Sqrt( lXiMomX*lXiMomX   + lXiMomY*lXiMomY );
1253         lXiTotMom       = TMath::Sqrt( lXiMomX*lXiMomX   + lXiMomY*lXiMomY   + lXiMomZ*lXiMomZ );
1254
1255                 lBachMomX = xi->MomBachX();
1256                 lBachMomY = xi->MomBachY();
1257                 lBachMomZ = xi->MomBachZ();
1258         lBachTransvMom  = TMath::Sqrt( lBachMomX*lBachMomX   + lBachMomY*lBachMomY );
1259         lBachTotMom     = TMath::Sqrt( lBachMomX*lBachMomX   + lBachMomY*lBachMomY  +  lBachMomZ*lBachMomZ  );
1260
1261
1262         lV0toXiCosineOfPointingAngle = xi->CosPointingAngle( xi->GetDecayVertexXi() );
1263
1264         lRapXi    = xi->RapXi();
1265         lRapOmega = xi->RapOmega();
1266         lEta      = xi->Eta();                          // Will not work ! need a method Pz(), Py() Px()
1267         lTheta    = xi->Theta() *180.0/TMath::Pi();     // in AODcascade.
1268         lPhi      = xi->Phi()   *180.0/TMath::Pi();     // Here, we will get eta, theta, phi for the V0 ...
1269         lAlphaXi  = xi->AlphaXi();
1270         lPtArmXi  = xi->PtArmXi();
1271
1272
1273                 // II.Step 7 - Azimuthal correlation study
1274                 //-------------
1275
1276         lTVect3MomXi.SetXYZ( lXiMomX, lXiMomY, lXiMomZ );
1277
1278         AliAODTrack *pTrackXi    = dynamic_cast<AliAODTrack*>( xi->GetDaughter(0) );
1279         AliAODTrack *nTrackXi    = dynamic_cast<AliAODTrack*>( xi->GetDaughter(1) );
1280         AliAODTrack *bachTrackXi = dynamic_cast<AliAODTrack*>( xi->GetDecayVertexXi()->GetDaughter(0) );
1281                 if (!pTrackXi || !nTrackXi || !bachTrackXi ) {
1282                         Printf("ERROR: Could not retrieve one of the 3 AOD daughter tracks of the cascade ...");
1283                         continue;
1284                 }
1285
1286         lArrTrackID[0] = pTrackXi   ->GetID();
1287         lArrTrackID[1] = nTrackXi   ->GetID();
1288         lArrTrackID[2] = bachTrackXi->GetID();
1289
1290   }// end of AOD treatment
1291
1292
1293   // -------------------------------------
1294   // III - Filling the TH1,2,3Fs
1295
1296
1297         // - III.Step 1
1298
1299         fHistVtxStatus                  ->Fill( lStatusTrackingPrimVtx   );  // 1 if tracking vtx = ok
1300
1301         if( lStatusTrackingPrimVtx ){
1302                 fHistPosTrkgPrimaryVtxX ->Fill( lTrkgPrimaryVtxPos[0]    );
1303                 fHistPosTrkgPrimaryVtxY ->Fill( lTrkgPrimaryVtxPos[1]    );
1304                 fHistPosTrkgPrimaryVtxZ ->Fill( lTrkgPrimaryVtxPos[2]    );
1305                 fHistTrkgPrimaryVtxRadius->Fill( lTrkgPrimaryVtxRadius3D );
1306         }
1307
1308         fHistPosBestPrimaryVtxX         ->Fill( lBestPrimaryVtxPos[0]    );
1309         fHistPosBestPrimaryVtxY         ->Fill( lBestPrimaryVtxPos[1]    );
1310         fHistPosBestPrimaryVtxZ         ->Fill( lBestPrimaryVtxPos[2]    );
1311         fHistBestPrimaryVtxRadius       ->Fill( lBestPrimaryVtxRadius3D  );
1312
1313         f2dHistTrkgPrimVtxVsBestPrimVtx->Fill( lTrkgPrimaryVtxRadius3D, lBestPrimaryVtxRadius3D );
1314
1315
1316         // - III.Step 2
1317         fHistEffMassXi                  ->Fill( lEffMassXi               );
1318         fHistChi2Xi                     ->Fill( lChi2Xi                  );     // Flag CascadeVtxer: Cut Variable a
1319         fHistDcaXiDaughters             ->Fill( lDcaXiDaughters          );     // Flag CascadeVtxer: Cut Variable e
1320         fHistDcaBachToPrimVertex        ->Fill( lDcaBachToPrimVertexXi   );     // Flag CascadeVtxer: Cut Variable d
1321         fHistXiCosineOfPointingAngle    ->Fill( lXiCosineOfPointingAngle );     // Flag CascadeVtxer: Cut Variable f
1322         fHistXiRadius                   ->Fill( lXiRadius                );     // Flag CascadeVtxer: Cut Variable g+h
1323
1324
1325         // - III.Step 3
1326         fHistMassLambdaAsCascDghter     ->Fill( lInvMassLambdaAsCascDghter );   // Flag CascadeVtxer: Cut Variable c
1327         fHistV0Chi2Xi                   ->Fill( lV0Chi2Xi                  );
1328         fHistDcaV0DaughtersXi           ->Fill( lDcaV0DaughtersXi          );
1329         fHistV0CosineOfPointingAngleXi  ->Fill( lV0CosineOfPointingAngleXi );
1330         fHistV0RadiusXi                 ->Fill( lV0RadiusXi                );
1331
1332         fHistDcaV0ToPrimVertexXi        ->Fill( lDcaV0ToPrimVertexXi       );   // Flag CascadeVtxer: Cut Variable b
1333         fHistDcaPosToPrimVertexXi       ->Fill( lDcaPosToPrimVertexXi      );
1334         fHistDcaNegToPrimVertexXi       ->Fill( lDcaNegToPrimVertexXi      );
1335
1336
1337         // - III.Step 4
1338         if( lChargeXi < 0 ){
1339                                         fHistMassXiMinus               ->Fill( lInvMassXiMinus    );
1340                                         fHistMassOmegaMinus            ->Fill( lInvMassOmegaMinus );
1341                 if(lIsBachelorPion)     fHistMassWithCombPIDXiMinus    ->Fill( lInvMassXiMinus    );
1342                 if(lIsBachelorKaon)     fHistMassWithCombPIDOmegaMinus ->Fill( lInvMassOmegaMinus );
1343         }
1344
1345         if( lChargeXi > 0 ){
1346                                         fHistMassXiPlus                ->Fill( lInvMassXiPlus     );
1347                                         fHistMassOmegaPlus             ->Fill( lInvMassOmegaPlus  );
1348                 if(lIsBachelorPion)     fHistMassWithCombPIDXiPlus     ->Fill( lInvMassXiPlus     );
1349                 if(lIsBachelorKaon)     fHistMassWithCombPIDOmegaPlus  ->Fill( lInvMassOmegaPlus  );
1350         }
1351
1352
1353         // - III.Step 6
1354         fHistXiTransvMom        ->Fill( lXiTransvMom   );
1355         fHistXiTotMom           ->Fill( lXiTotMom      );
1356
1357         fHistBachTransvMom      ->Fill( lBachTransvMom );
1358         fHistBachTotMom         ->Fill( lBachTotMom    );
1359
1360         fHistChargeXi           ->Fill( lChargeXi      );
1361         fHistV0toXiCosineOfPointingAngle->Fill( lV0toXiCosineOfPointingAngle );
1362
1363         fHistRapXi              ->Fill( lRapXi               );
1364         fHistRapOmega           ->Fill( lRapOmega            );
1365         fHistEta                ->Fill( lEta                 );
1366         fHistTheta              ->Fill( lTheta               );
1367         fHistPhi                ->Fill( lPhi                 );
1368
1369         f2dHistArmenteros       ->Fill( lAlphaXi, lPtArmXi   );
1370
1371         // - III.Step 5
1372         if( lChargeXi < 0 ) {
1373                 f2dHistEffMassLambdaVsEffMassXiMinus->Fill( lInvMassLambdaAsCascDghter, lInvMassXiMinus );
1374                 f2dHistEffMassXiVsEffMassOmegaMinus ->Fill( lInvMassXiMinus, lInvMassOmegaMinus );
1375                 f2dHistXiRadiusVsEffMassXiMinus     ->Fill( lXiRadius, lInvMassXiMinus );
1376                 f2dHistXiRadiusVsEffMassOmegaMinus  ->Fill( lXiRadius, lInvMassOmegaMinus );
1377                 f3dHistXiPtVsEffMassVsYXiMinus      ->Fill( lXiTransvMom, lInvMassXiMinus,    lRapXi    );
1378                 f3dHistXiPtVsEffMassVsYOmegaMinus   ->Fill( lXiTransvMom, lInvMassOmegaMinus, lRapOmega );
1379                 if(lIsPosInXiProton)                            f3dHistXiPtVsEffMassVsYWithCombPIDXiMinus     ->Fill( lXiTransvMom, lInvMassXiMinus,    lRapXi    );
1380                 if(lIsBachelorKaon)                             f3dHistXiPtVsEffMassVsYWithCombPIDOmegaMinus  ->Fill( lXiTransvMom, lInvMassOmegaMinus, lRapOmega );
1381                 if(lIsBachelorKaon && lIsPosInOmegaProton)      f3dHistXiPtVsEffMassVsYWith2CombPIDOmegaMinus ->Fill( lXiTransvMom, lInvMassOmegaMinus, lRapOmega );
1382         }
1383         else{
1384                 f2dHistEffMassLambdaVsEffMassXiPlus ->Fill( lInvMassLambdaAsCascDghter, lInvMassXiPlus );
1385                 f2dHistEffMassXiVsEffMassOmegaPlus  ->Fill( lInvMassXiPlus, lInvMassOmegaPlus );
1386                 f2dHistXiRadiusVsEffMassXiPlus      ->Fill( lXiRadius, lInvMassXiPlus);
1387                 f2dHistXiRadiusVsEffMassOmegaPlus   ->Fill( lXiRadius, lInvMassOmegaPlus );
1388                 f3dHistXiPtVsEffMassVsYXiPlus       ->Fill( lXiTransvMom, lInvMassXiPlus,    lRapXi    );
1389                 f3dHistXiPtVsEffMassVsYOmegaPlus    ->Fill( lXiTransvMom, lInvMassOmegaPlus, lRapOmega );
1390                 if(lIsNegInXiProton)                            f3dHistXiPtVsEffMassVsYWithCombPIDXiPlus     ->Fill( lXiTransvMom, lInvMassXiPlus,    lRapXi    );
1391                 if(lIsBachelorKaon )                            f3dHistXiPtVsEffMassVsYWithCombPIDOmegaPlus  ->Fill( lXiTransvMom, lInvMassOmegaPlus, lRapOmega );
1392                 if(lIsBachelorKaon && lIsNegInOmegaProton)      f3dHistXiPtVsEffMassVsYWith2CombPIDOmegaPlus ->Fill( lXiTransvMom, lInvMassOmegaPlus, lRapOmega );
1393         }
1394
1395
1396
1397
1398
1399         // - III.Step 7
1400
1401         if( lChargeXi < 0 ){
1402                 DoAngularCorrelation("Xi-",    lInvMassXiMinus,    lArrTrackID, lTVect3MomXi, lEta);
1403                 DoAngularCorrelation("Omega-", lInvMassOmegaMinus, lArrTrackID, lTVect3MomXi, lEta);
1404         }
1405         else{
1406                 DoAngularCorrelation("Xi+",    lInvMassXiPlus,    lArrTrackID, lTVect3MomXi, lEta);
1407                 DoAngularCorrelation("Omega+", lInvMassOmegaPlus, lArrTrackID, lTVect3MomXi, lEta);
1408         }
1409
1410
1411     }// end of the Cascade loop (ESD or AOD)
1412
1413
1414   // Post output data.
1415  PostData(1, fListHistCascade);
1416 }
1417
1418
1419 void AliAnalysisTaskCheckCascade::DoAngularCorrelation( const Char_t   *lCascType,
1420                                                               Double_t  lInvMassCascade,
1421                                                         const Int_t    *lArrTrackID,
1422                                                               TVector3 &lTVect3MomXi,
1423                                                               Double_t  lEtaXi ){
1424   // Perform the Delta(Phi)Delta(Eta) analysis
1425   // by properly filling the THnSparseF
1426
1427         TString lStrCascType( lCascType );
1428
1429         Double_t lCascPdgMass = 0.0;
1430         if( lStrCascType.Contains("Xi") )       lCascPdgMass = 1.3217;
1431         if( lStrCascType.Contains("Omega") )    lCascPdgMass = 1.6724;
1432
1433         if( lInvMassCascade > lCascPdgMass + 0.010) return;
1434         if( lInvMassCascade < lCascPdgMass - 0.010) return;
1435         // Check the Xi- candidate is within the proper mass window m0 +- 10 MeV
1436
1437
1438         // 1st loop: check there is no track with a higher pt ...
1439         // = The cascade is meant to be a leading particle : Pt(Casc) > any track in the event
1440         for(Int_t TrckIdx = 0; TrckIdx < (InputEvent())->GetNumberOfTracks() ; TrckIdx++ )
1441         {// Loop over all the tracks of the event
1442
1443                 AliVTrack *lCurrentTrck = dynamic_cast<AliVTrack*>( (InputEvent())->GetTrack( TrckIdx ) );
1444                         if (!lCurrentTrck ) {
1445                                 Printf("ERROR Correl. Study : Could not retrieve a track while looping over the event tracks ...");
1446                                 continue;
1447                         }
1448                 if(lTVect3MomXi.Pt() < lCurrentTrck->Pt() ) return;
1449                 // Room for improvement: //FIXME
1450                 // 1. There is a given resolution on pt : maybe release the cut Pt(casc) < Pt(track)*90% ?
1451                 // 2. Apply this cut only when DeltaPhi(casc, track) > 90 deg = when track is in the away-side ?
1452
1453         }// end control loop
1454
1455         // 2nd loop: filling loop
1456         for(Int_t TrckIdx = 0; TrckIdx < (InputEvent())->GetNumberOfTracks() ; TrckIdx++ )
1457         {// Loop over all the tracks of the event
1458
1459                 AliVTrack *lCurrentTrck = dynamic_cast<AliVTrack*>( (InputEvent())->GetTrack( TrckIdx ) );
1460                         if (!lCurrentTrck ) {
1461                                 Printf("ERROR Correl. Study : Could not retrieve a track while looping over the event tracks ...");
1462                                 continue;
1463                         }
1464
1465                 // Room for improvement: //FIXME
1466                 // 1. Loop only on primary tracks ?
1467                 // 2. Exclude the tracks that build the condisdered cascade = the bachelor + the V0 dghters
1468                 //     This may bias the outcome, especially for low multplicity events.
1469                 // Note : For ESD event, track ID == track index.
1470                         if(lCurrentTrck->GetID() == lArrTrackID[0]) continue;
1471                         if(lCurrentTrck->GetID() == lArrTrackID[1]) continue;
1472                         if(lCurrentTrck->GetID() == lArrTrackID[2]) continue;
1473
1474                 TVector3 lTVect3MomTrck(lCurrentTrck->Px(), lCurrentTrck->Py(), lCurrentTrck->Pz() );
1475
1476                 // 2 hypotheses made here :
1477                 //   - The Xi trajectory is a straight line,
1478                 //   - The Xi doesn't loose any energy by crossing the first layer(s) of ITS, if ever;
1479                 //      So, meaning hyp: vect p(Xi) at the emission = vect p(Xi) at the decay vertex
1480                 //      By doing this, we introduce a systematic error on the cascade Phi ...
1481                 // Room for improvement: take into account the curvature of the Xi trajectory //FIXME
1482
1483                 Double_t lHnSpFillVar[5] = {0.};
1484                 lHnSpFillVar[0] = lTVect3MomXi.DeltaPhi(lTVect3MomTrck) * 180.0/TMath::Pi(); // Delta phi(Casc,Track) (deg)
1485                 if(lHnSpFillVar[0] < -50.0) lHnSpFillVar[0] += 360.0;
1486                 lHnSpFillVar[1] = lEtaXi - lCurrentTrck->Eta();                            // Delta eta(Casc,Track)
1487                 lHnSpFillVar[2] = lTVect3MomXi.Pt();                                       // Pt_{Casc}
1488                 lHnSpFillVar[3] = lCurrentTrck->Pt();                                      // Pt_{any track}
1489                 lHnSpFillVar[4] = lInvMassCascade;                                         // Eff. Inv Mass (control var)
1490
1491                 if(      lStrCascType.Contains("Xi-") )      fHnSpAngularCorrXiMinus    ->Fill( lHnSpFillVar );
1492                 else if( lStrCascType.Contains("Xi+") )      fHnSpAngularCorrXiPlus     ->Fill( lHnSpFillVar );
1493                 else if( lStrCascType.Contains("Omega-") )   fHnSpAngularCorrOmegaMinus ->Fill( lHnSpFillVar );
1494                 else if( lStrCascType.Contains("Omega+") )   fHnSpAngularCorrOmegaPlus  ->Fill( lHnSpFillVar );
1495
1496         }// end - Loop over all the tracks in the event
1497
1498 }
1499
1500
1501
1502
1503
1504
1505
1506 //________________________________________________________________________
1507 void AliAnalysisTaskCheckCascade::Terminate(Option_t *)
1508 {
1509   // Draw result to the screen
1510   // Called once at the end of the query
1511
1512   TList *cRetrievedList = 0x0;
1513          cRetrievedList = (TList*)GetOutputData(1);
1514         if(!cRetrievedList){
1515                 Printf("ERROR - AliAnalysisTaskCheckCascade: ouput data container list not available\n");
1516                 return;
1517         }
1518
1519   fHistTrackMultiplicity = dynamic_cast<TH1F*> (   cRetrievedList->FindObject("fHistTrackMultiplicity") );
1520         if (!fHistTrackMultiplicity) {
1521                 Printf("ERROR - AliAnalysisTaskCheckCascade: fHistTrackMultiplicity not available\n");
1522                 return;
1523         }
1524
1525   fHistCascadeMultiplicity = dynamic_cast<TH1F*> ( cRetrievedList->FindObject("fHistCascadeMultiplicity"));
1526         if (!fHistCascadeMultiplicity) {
1527                 Printf("ERROR - AliAnalysisTaskCheckCascade: fHistCascadeMultiplicity not available\n");
1528                 return;
1529         }
1530
1531   TCanvas *canCheckCascade = new TCanvas("AliAnalysisTaskCheckCascade","Multiplicity",10,10,510,510);
1532   canCheckCascade->cd(1)->SetLogy();
1533
1534   fHistTrackMultiplicity->SetMarkerStyle(22);
1535   fHistTrackMultiplicity->DrawCopy("E");
1536
1537   fHistCascadeMultiplicity->SetMarkerStyle(26);
1538   fHistCascadeMultiplicity->DrawCopy("ESAME");
1539
1540 }