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 //            This task is for QAing the Cascades from ESD and AOD
  18 //            Origin   : Antonin Maire Fev2008, antonin.maire@ires.in2p3.fr
  19 //            Modified : A.Maire June 2009
  20 //-----------------------------------------------------------------
  21
  22
  23
  24
  25 class TTree;
  26 class TParticle;
  27 class TVector3;
  28
  29 //class AliMCEventHandler;
  30 //class AliMCEvent;
  31 //class AliStack;
  32
  33 class AliESDVertex;
  34 class AliAODVertex;
  35 class AliESDv0;
  36 class AliAODv0;
  37
  38 #include <Riostream.h>
  39
  40 #include "TList.h"
  41 #include "TH1.h"
  42 #include "TH2.h"
  43 #include "TCanvas.h"
  44 #include "TMath.h"
  45
  46 #include "AliLog.h"
  47
  48 #include "AliESDEvent.h"
  49 #include "AliAODEvent.h"
  50 //#include "AliCascadeVertexer.h"
  51
  52 #include "AliESDcascade.h"
  53 #include "AliAODcascade.h"
  54
  55 #include "AliAnalysisTaskCheckCascade.h"
  56
  57 ClassImp(AliAnalysisTaskCheckCascade)
  58
  59
  60
  61 //________________________________________________________________________
  62 AliAnalysisTaskCheckCascade::AliAnalysisTaskCheckCascade()
  63   : AliAnalysisTaskSE(), fAnalysisType("ESD"), fCollidingSystems(0),
  64
  65         // - Cascade part initialisation
  66     fListHistCascade(0),
  67     fHistTrackMultiplicity(0), fHistCascadeMultiplicity(0),
  68     fHistVtxStatus(0),
  69
  70     fHistPosTrkgPrimaryVtxX(0), fHistPosTrkgPrimaryVtxY(0), fHistPosTrkgPrimaryVtxZ(0), fHistTrkgPrimaryVtxRadius(0),
  71     fHistPosBestPrimaryVtxX(0), fHistPosBestPrimaryVtxY(0), fHistPosBestPrimaryVtxZ(0), fHistBestPrimaryVtxRadius(0),
  72     f2dHistTrkgPrimVtxVsBestPrimVtx(0),
  73
  74     fHistEffMassXi(0),  fHistChi2Xi(0),
  75     fHistDcaXiDaughters(0), fHistDcaBachToPrimVertex(0), fHistXiCosineOfPointingAngle(0), fHistXiRadius(0),
  76
  77     fHistMassLambdaAsCascDghter(0),
  78     fHistV0Chi2Xi(0),
  79     fHistDcaV0DaughtersXi(0),
  80     fHistDcaV0ToPrimVertexXi(0),
  81     fHistV0CosineOfPointingAngleXi(0),
  82     fHistV0RadiusXi(0),
  83     fHistDcaPosToPrimVertexXi(0), fHistDcaNegToPrimVertexXi(0),
  84
  85     fHistMassXiMinus(0), fHistMassXiPlus(0),
  86     fHistMassOmegaMinus(0), fHistMassOmegaPlus(0),
  87     fHistMassWithCombPIDXiMinus(0), fHistMassWithCombPIDXiPlus(0),
  88     fHistMassWithCombPIDOmegaMinus(0), fHistMassWithCombPIDOmegaPlus(0),
  89
  90     fHistXiTransvMom(0),    fHistXiTotMom(0),
  91     fHistBachTransvMom(0),   fHistBachTotMom(0),
  92
  93     fHistChargeXi(0),
  94     fHistV0toXiCosineOfPointingAngle(0),
  95
  96     fHistRapXi(0), fHistRapOmega(0), fHistEta(0),
  97     fHistTheta(0), fHistPhi(0),
  98
  99     f2dHistArmenteros(0),
 100     f2dHistEffMassLambdaVsEffMassXiMinus(0), f2dHistEffMassXiVsEffMassOmegaMinus(0),
 101     f2dHistEffMassLambdaVsEffMassXiPlus(0), f2dHistEffMassXiVsEffMassOmegaPlus(0),
 102     f2dHistXiRadiusVsEffMassXiMinus(0), f2dHistXiRadiusVsEffMassXiPlus(0),
 103     f2dHistXiRadiusVsEffMassOmegaMinus(0), f2dHistXiRadiusVsEffMassOmegaPlus(0),
 104
 105     f2dHistXiPtVsEffMassXiMinus(0), f2dHistXiPtVsEffMassXiPlus(0),
 106     f2dHistXiPtVsEffMassOmegaMinus(0), f2dHistXiPtVsEffMassOmegaPlus(0)
 107
 108 {
 109   // Dummy Constructor
 110 }
 111
 112
 113
 114
 115
 116
 117
 118
 119 //________________________________________________________________________
 120 AliAnalysisTaskCheckCascade::AliAnalysisTaskCheckCascade(const char *name)
 121   : AliAnalysisTaskSE(name), fAnalysisType("ESD"), fCollidingSystems(0),
 122
 123         // - Cascade part initialisation
 124     fListHistCascade(0),
 125     fHistTrackMultiplicity(0), fHistCascadeMultiplicity(0),
 126     fHistVtxStatus(0),
 127
 128     fHistPosTrkgPrimaryVtxX(0), fHistPosTrkgPrimaryVtxY(0), fHistPosTrkgPrimaryVtxZ(0), fHistTrkgPrimaryVtxRadius(0),
 129     fHistPosBestPrimaryVtxX(0), fHistPosBestPrimaryVtxY(0), fHistPosBestPrimaryVtxZ(0), fHistBestPrimaryVtxRadius(0),
 130     f2dHistTrkgPrimVtxVsBestPrimVtx(0),
 131
 132     fHistEffMassXi(0),  fHistChi2Xi(0),
 133     fHistDcaXiDaughters(0), fHistDcaBachToPrimVertex(0), fHistXiCosineOfPointingAngle(0), fHistXiRadius(0),
 134
 135     fHistMassLambdaAsCascDghter(0),
 136     fHistV0Chi2Xi(0),
 137     fHistDcaV0DaughtersXi(0),
 138     fHistDcaV0ToPrimVertexXi(0),
 139     fHistV0CosineOfPointingAngleXi(0),
 140     fHistV0RadiusXi(0),
 141     fHistDcaPosToPrimVertexXi(0), fHistDcaNegToPrimVertexXi(0),
 142
 143     fHistMassXiMinus(0), fHistMassXiPlus(0),
 144     fHistMassOmegaMinus(0), fHistMassOmegaPlus(0),
 145     fHistMassWithCombPIDXiMinus(0), fHistMassWithCombPIDXiPlus(0),
 146     fHistMassWithCombPIDOmegaMinus(0), fHistMassWithCombPIDOmegaPlus(0),
 147
 148     fHistXiTransvMom(0),    fHistXiTotMom(0),
 149     fHistBachTransvMom(0),   fHistBachTotMom(0),
 150
 151     fHistChargeXi(0),
 152     fHistV0toXiCosineOfPointingAngle(0),
 153
 154     fHistRapXi(0), fHistRapOmega(0), fHistEta(0),
 155     fHistTheta(0), fHistPhi(0),
 156
 157     f2dHistArmenteros(0),
 158     f2dHistEffMassLambdaVsEffMassXiMinus(0), f2dHistEffMassXiVsEffMassOmegaMinus(0),
 159     f2dHistEffMassLambdaVsEffMassXiPlus(0), f2dHistEffMassXiVsEffMassOmegaPlus(0),
 160     f2dHistXiRadiusVsEffMassXiMinus(0), f2dHistXiRadiusVsEffMassXiPlus(0),
 161     f2dHistXiRadiusVsEffMassOmegaMinus(0), f2dHistXiRadiusVsEffMassOmegaPlus(0),
 162
 163     f2dHistXiPtVsEffMassXiMinus(0), f2dHistXiPtVsEffMassXiPlus(0),
 164     f2dHistXiPtVsEffMassOmegaMinus(0), f2dHistXiPtVsEffMassOmegaPlus(0)
 165
 166 {
 167   // Constructor
 168
 169   // Define input and output slots here
 170   // Input slot #0 works with a TChain
 171
 172   // Output slot #0 writes into a TList container (Cascade)
 173   DefineOutput(1, TList::Class());
 174 }
 175
 176
 177
 178
 179
 180
 181 //________________________________________________________________________
 182 void AliAnalysisTaskCheckCascade::UserCreateOutputObjects()
 183 {
 184   // Create histograms
 185   // Called once
 186
 187
 188
 189  fListHistCascade = new TList();
 190
 191
 192         // - General histos
 193
 194 if(! fHistTrackMultiplicity) {
 195         if(fCollidingSystems)// AA collisions
 196                 fHistTrackMultiplicity = new TH1F("fHistTrackMultiplicity",
 197                         "Multiplicity distribution;Number of tracks;Events",
 198                         200, 0, 40000);
 199         else // pp collisions
 200                 fHistTrackMultiplicity = new TH1F("fHistTrackMultiplicity",
 201                         "Track Multiplicity;Nbr of tracks/Evt;Events",
 202                         200, 0, 200);
 203         fListHistCascade->Add(fHistTrackMultiplicity);
 204 }
 205
 206 if(! fHistCascadeMultiplicity) {
 207         if(fCollidingSystems)// AA collisions
 208                 fHistCascadeMultiplicity = new TH1F("fHistCascadeMultiplicity",
 209                         "Multiplicity distribution;Number of Cascades;Events",
 210                         25, 0, 25);
 211         else // pp collisions
 212                 fHistCascadeMultiplicity = new TH1F("fHistCascadeMultiplicity",
 213                         "Cascades per event;Nbr of Cascades/Evt;Events",
 214                         10, 0, 10);
 215         fListHistCascade->Add(fHistCascadeMultiplicity);
 216 }
 217
 218
 219
 220 if(! fHistVtxStatus ){
 221         fHistVtxStatus   = new TH1F( "fHistVtxStatus" , "Does a Trckg Prim.vtx exist ?; true=1 or false=0; Nb of Events" , 4, -1.0, 3.0 );
 222         fListHistCascade->Add(fHistVtxStatus);
 223 }
 224
 225
 226
 227
 228
 229
 230         // - Vertex Positions
 231
 232 if(! fHistPosTrkgPrimaryVtxX ){
 233         fHistPosTrkgPrimaryVtxX   = new TH1F( "fHistPosTrkgPrimaryVtxX" , "Trkg Prim. Vertex Position in x; x (cm); Events" , 200, -0.5, 0.5 );
 234         fListHistCascade->Add(fHistPosTrkgPrimaryVtxX);
 235 }
 236
 237
 238 if(! fHistPosTrkgPrimaryVtxY){
 239         fHistPosTrkgPrimaryVtxY   = new TH1F( "fHistPosTrkgPrimaryVtxY" , "Trkg Prim. Vertex Position in y; y (cm); Events" , 200, -0.5, 0.5 );
 240         fListHistCascade->Add(fHistPosTrkgPrimaryVtxY);
 241 }
 242
 243 if(! fHistPosTrkgPrimaryVtxZ ){
 244         fHistPosTrkgPrimaryVtxZ   = new TH1F( "fHistPosTrkgPrimaryVtxZ" , "Trkg Prim. Vertex Position in z; z (cm); Events" , 100, -15.0, 15.0 );
 245         fListHistCascade->Add(fHistPosTrkgPrimaryVtxZ);
 246 }
 247
 248 if(! fHistTrkgPrimaryVtxRadius ){
 249         fHistTrkgPrimaryVtxRadius  = new TH1F( "fHistTrkgPrimaryVtxRadius",  "Trkg Prim. Vertex radius; r (cm); Events" , 150, 0., 15.0 );
 250         fListHistCascade->Add(fHistTrkgPrimaryVtxRadius);
 251 }
 252
 253
 254
 255
 256 if(! fHistPosBestPrimaryVtxX ){
 257         fHistPosBestPrimaryVtxX   = new TH1F( "fHistPosBestPrimaryVtxX" , "Best Prim. Vertex Position in x; x (cm); Events" , 200, -0.5, 0.5 );
 258         fListHistCascade->Add(fHistPosBestPrimaryVtxX);
 259 }
 260
 261 if(! fHistPosBestPrimaryVtxY){
 262         fHistPosBestPrimaryVtxY   = new TH1F( "fHistPosBestPrimaryVtxY" , "Best Prim. Vertex Position in y; y (cm); Events" , 200, -0.5, 0.5 );
 263         fListHistCascade->Add(fHistPosBestPrimaryVtxY);
 264 }
 265
 266 if(! fHistPosBestPrimaryVtxZ ){
 267         fHistPosBestPrimaryVtxZ   = new TH1F( "fHistPosBestPrimaryVtxZ" , "Best Prim. Vertex Position in z; z (cm); Events" , 100, -15.0, 15.0 );
 268         fListHistCascade->Add(fHistPosBestPrimaryVtxZ);
 269 }
 270
 271 if(! fHistBestPrimaryVtxRadius ){
 272         fHistBestPrimaryVtxRadius  = new TH1F( "fHistBestPrimaryVtxRadius",  "Best Prim.  vertex radius; r (cm); Events" , 150, 0., 15.0 );
 273         fListHistCascade->Add(fHistBestPrimaryVtxRadius);
 274 }
 275
 276 if(! f2dHistTrkgPrimVtxVsBestPrimVtx) {
 277         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.);
 278         fListHistCascade->Add(f2dHistTrkgPrimVtxVsBestPrimVtx);
 279 }
 280
 281
 282
 283
 284 // - Typical histos for cascades
 285
 286
 287 if(! fHistEffMassXi) {
 288      fHistEffMassXi = new TH1F("fHistEffMassXi", "Cascade candidates ; Invariant Mass (GeV/c^{2}) ; Counts", 200, 1.2, 2.0);
 289      fListHistCascade->Add(fHistEffMassXi);
 290 }
 291
 292 if(! fHistChi2Xi ){
 293         fHistChi2Xi = new TH1F("fHistChi2Xi", "Cascade #chi^{2}; #chi^{2}; Number of Cascades", 160, 0, 40);
 294         fListHistCascade->Add(fHistChi2Xi);
 295 }
 296
 297 if(! fHistDcaXiDaughters ){
 298         fHistDcaXiDaughters = new TH1F( "fHistDcaXiDaughters",  "DCA between Xi Daughters; DCA (cm) ; Number of Cascades", 100, 0., 0.5);
 299         fListHistCascade->Add(fHistDcaXiDaughters);
 300 }
 301
 302 if(! fHistDcaBachToPrimVertex) {
 303         fHistDcaBachToPrimVertex = new TH1F("fHistDcaBachToPrimVertex", "DCA of Bach. to Prim. Vertex;DCA (cm);Number of Cascades", 250, 0., 0.25);
 304         fListHistCascade->Add(fHistDcaBachToPrimVertex);
 305 }
 306
 307 if(! fHistXiCosineOfPointingAngle) {
 308         fHistXiCosineOfPointingAngle = new TH1F("fHistXiCosineOfPointingAngle", "Cosine of Xi Pointing Angle; Cos (Xi Point.Angl);Number of Xis", 200, 0.98, 1.0);
 309         fListHistCascade->Add(fHistXiCosineOfPointingAngle);
 310 }
 311
 312 if(! fHistXiRadius ){
 313         fHistXiRadius  = new TH1F( "fHistXiRadius",  "Casc. decay transv. radius; r (cm); Counts" , 200, 0., 20.0 );
 314         fListHistCascade->Add(fHistXiRadius);
 315 }
 316
 317
 318 // - Histos about ~ the "V0 part" of the cascade,  coming by inheritance from AliESDv0
 319
 320
 321
 322 if (! fHistMassLambdaAsCascDghter) {
 323      fHistMassLambdaAsCascDghter = new TH1F("fHistMassLambdaAsCascDghter","#Lambda associated to Casc. candidates;Eff. Mass (GeV/c^{2});Counts", 160,1.00,1.8);
 324     fListHistCascade->Add(fHistMassLambdaAsCascDghter);
 325 }
 326
 327 if (! fHistV0Chi2Xi) {
 328         fHistV0Chi2Xi = new TH1F("fHistV0Chi2Xi", "V0 #chi^{2}, in cascade; #chi^{2};Counts", 160, 0, 40);
 329         fListHistCascade->Add(fHistV0Chi2Xi);
 330 }
 331
 332 if (! fHistDcaV0DaughtersXi) {
 333         fHistDcaV0DaughtersXi = new TH1F("fHistDcaV0DaughtersXi", "DCA between V0 daughters, in cascade;DCA (cm);Number of V0s", 120, 0., 0.6);
 334         fListHistCascade->Add(fHistDcaV0DaughtersXi);
 335 }
 336
 337 if (! fHistDcaV0ToPrimVertexXi) {
 338         fHistDcaV0ToPrimVertexXi = new TH1F("fHistDcaV0ToPrimVertexXi", "DCA of V0 to Prim. Vertex, in cascade;DCA (cm);Number of Cascades", 200, 0., 1.);
 339         fListHistCascade->Add(fHistDcaV0ToPrimVertexXi);
 340 }
 341
 342 if (! fHistV0CosineOfPointingAngleXi) {
 343         fHistV0CosineOfPointingAngleXi = new TH1F("fHistV0CosineOfPointingAngleXi", "Cosine of V0 Pointing Angle, in cascade;Cos(V0 Point. Angl); Counts", 200, 0.98, 1.0);
 344         fListHistCascade->Add(fHistV0CosineOfPointingAngleXi);
 345 }
 346
 347 if (! fHistV0RadiusXi) {
 348         fHistV0RadiusXi = new TH1F("fHistV0RadiusXi", "V0 decay radius, in cascade; radius (cm); Counts", 200, 0, 20);
 349         fListHistCascade->Add(fHistV0RadiusXi);
 350 }
 351
 352 if (! fHistDcaPosToPrimVertexXi) {
 353         fHistDcaPosToPrimVertexXi = new TH1F("fHistDcaPosToPrimVertexXi", "DCA of V0 pos daughter to Prim. Vertex;DCA (cm);Counts", 300, 0, 3);
 354         fListHistCascade->Add(fHistDcaPosToPrimVertexXi);
 355 }
 356
 357 if (! fHistDcaNegToPrimVertexXi) {
 358         fHistDcaNegToPrimVertexXi = new TH1F("fHistDcaNegToPrimVertexXi", "DCA of V0 neg daughter to Prim. Vertex;DCA (cm);Counts", 300, 0, 3);
 359         fListHistCascade->Add(fHistDcaNegToPrimVertexXi);
 360 }
 361
 362
 363
 364
 365         // - Effective mass histos for cascades.
 366 // By cascade hyp
 367 if (! fHistMassXiMinus) {
 368     fHistMassXiMinus = new TH1F("fHistMassXiMinus","#Xi^{-} candidates;M( #Lambda , #pi^{-} ) (GeV/c^{2});Counts", 200,1.2,2.0);
 369     fListHistCascade->Add(fHistMassXiMinus);
 370 }
 371
 372 if (! fHistMassXiPlus) {
 373     fHistMassXiPlus = new TH1F("fHistMassXiPlus","#Xi^{+} candidates;M( #bar{#Lambda}^{0} , #pi^{+} ) (GeV/c^{2});Counts",200,1.2,2.0);
 374     fListHistCascade->Add(fHistMassXiPlus);
 375 }
 376
 377 if (! fHistMassOmegaMinus) {
 378         fHistMassOmegaMinus = new TH1F("fHistMassOmegaMinus","#Omega^{-} candidates;M( #Lambda , K^{-} ) (GeV/c^{2});Counts", 250,1.5,2.5);
 379     fListHistCascade->Add(fHistMassOmegaMinus);
 380 }
 381
 382 if (! fHistMassOmegaPlus) {
 383         fHistMassOmegaPlus = new TH1F("fHistMassOmegaPlus","#Omega^{+} candidates;M( #bar{#Lambda}^{0} , K^{+} ) (GeV/c^{2});Counts", 250,1.5,2.5);
 384     fListHistCascade->Add(fHistMassOmegaPlus);
 385 }
 386
 387 // By cascade hyp + bachelor PID
 388 if (! fHistMassWithCombPIDXiMinus) {
 389     fHistMassWithCombPIDXiMinus = new TH1F("fHistMassWithCombPIDXiMinus","#Xi^{-} candidates, with Bach. comb. PID;M( #Lambda , #pi^{-} ) (GeV/c^{2});Counts", 200,1.2,2.0);
 390     fListHistCascade->Add(fHistMassWithCombPIDXiMinus);
 391 }
 392
 393 if (! fHistMassWithCombPIDXiPlus) {
 394     fHistMassWithCombPIDXiPlus = new TH1F("fHistMassWithCombPIDXiPlus","#Xi^{+} candidates, with Bach. comb. PID;M( #bar{#Lambda}^{0} , #pi^{+} ) (GeV/c^{2});Counts",200,1.2,2.0);
 395     fListHistCascade->Add(fHistMassWithCombPIDXiPlus);
 396 }
 397
 398 if (! fHistMassWithCombPIDOmegaMinus) {
 399         fHistMassWithCombPIDOmegaMinus = new TH1F("fHistMassWithCombPIDOmegaMinus","#Omega^{-} candidates, with Bach. comb. PID;M( #Lambda , K^{-} ) (GeV/c^{2});Counts", 250,1.5,2.5);
 400     fListHistCascade->Add(fHistMassWithCombPIDOmegaMinus);
 401 }
 402
 403 if (! fHistMassWithCombPIDOmegaPlus) {
 404         fHistMassWithCombPIDOmegaPlus = new TH1F("fHistMassWithCombPIDOmegaPlus","#Omega^{+} candidates, with Bach. comb. PID;M( #bar{#Lambda}^{0} , K^{+} ) (GeV/c^{2});Counts", 250,1.5,2.5);
 405     fListHistCascade->Add(fHistMassWithCombPIDOmegaPlus);
 406 }
 407
 408
 409
 410         // - Complements for QA
 411
 412 if(! fHistXiTransvMom ){
 413         fHistXiTransvMom  = new TH1F( "fHistXiTransvMom" , "Xi transverse momentum ; p_{t}(#Xi) (GeV/c); Counts", 100, 0.0, 10.0);
 414         fListHistCascade->Add(fHistXiTransvMom);
 415 }
 416
 417 if(! fHistXiTotMom ){
 418         fHistXiTotMom  = new TH1F( "fHistXiTotMom" , "Xi momentum norm; p_{tot}(#Xi) (GeV/c); Counts", 150, 0.0, 15.0);
 419         fListHistCascade->Add(fHistXiTotMom);
 420 }
 421
 422
 423 if(! fHistBachTransvMom ){
 424         fHistBachTransvMom  = new TH1F( "fHistBachTransvMom" , "Bach. transverse momentum ; p_{t}(Bach.) (GeV/c); Counts", 100, 0.0, 5.0);
 425         fListHistCascade->Add(fHistBachTransvMom);
 426 }
 427
 428 if(! fHistBachTotMom ){
 429         fHistBachTotMom  = new TH1F( "fHistBachTotMom" , "Bach. momentum norm; p_{tot}(Bach.) (GeV/c); Counts", 100, 0.0, 5.0);
 430         fListHistCascade->Add(fHistBachTotMom);
 431 }
 432
 433
 434 if(! fHistChargeXi ){
 435         fHistChargeXi  = new TH1F( "fHistChargeXi" , "Charge of casc. candidates ; Sign ; Counts", 5, -2.0, 3.0);
 436         fListHistCascade->Add(fHistChargeXi);
 437 }
 438
 439
 440 if (! fHistV0toXiCosineOfPointingAngle) {
 441         fHistV0toXiCosineOfPointingAngle = new TH1F("fHistV0toXiCosineOfPointingAngle", "Cos. of V0 Ptng Angl / Xi vtx ;Cos(V0 Point. Angl / Xi vtx); Counts", 100, 0.99, 1.0);
 442         fListHistCascade->Add(fHistV0toXiCosineOfPointingAngle);
 443 }
 444
 445
 446 if(! fHistRapXi ){
 447         fHistRapXi  = new TH1F( "fHistRapXi" , "Rapidity of Xi candidates ; y ; Counts", 200, -5.0, 5.0);
 448         fListHistCascade->Add(fHistRapXi);
 449 }
 450
 451 if(! fHistRapOmega ){
 452         fHistRapOmega  = new TH1F( "fHistRapOmega" , "Rapidity of Omega candidates ; y ; Counts", 200, -5.0, 5.0);
 453         fListHistCascade->Add(fHistRapOmega);
 454 }
 455
 456 if(! fHistEta ){
 457         fHistEta  = new TH1F( "fHistEta" , "Pseudo-rap. of casc. candidates ; #eta ; Counts", 120, -3.0, 3.0);
 458         fListHistCascade->Add(fHistEta);
 459 }
 460
 461 if(! fHistTheta ){
 462         fHistTheta  = new TH1F( "fHistTheta" , "#theta of casc. candidates ; #theta (deg) ; Counts", 180, 0., 180.0);
 463         fListHistCascade->Add(fHistTheta);
 464 }
 465
 466 if(! fHistPhi ){
 467         fHistPhi  = new TH1F( "fHistPhi" , "#phi of casc. candidates ; #phi (deg) ; Counts", 360, 0., 360.);
 468         fListHistCascade->Add(fHistPhi);
 469 }
 470
 471
 472 if(! f2dHistArmenteros) {
 473         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);
 474         fListHistCascade->Add(f2dHistArmenteros);
 475 }
 476
 477 //-------
 478
 479 if(! f2dHistEffMassLambdaVsEffMassXiMinus) {
 480         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);
 481         fListHistCascade->Add(f2dHistEffMassLambdaVsEffMassXiMinus);
 482 }
 483
 484 if(! f2dHistEffMassXiVsEffMassOmegaMinus) {
 485         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);
 486         fListHistCascade->Add(f2dHistEffMassXiVsEffMassOmegaMinus);
 487 }
 488
 489 if(! f2dHistEffMassLambdaVsEffMassXiPlus) {
 490         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);
 491         fListHistCascade->Add(f2dHistEffMassLambdaVsEffMassXiPlus);
 492 }
 493
 494 if(! f2dHistEffMassXiVsEffMassOmegaPlus) {
 495         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);
 496         fListHistCascade->Add(f2dHistEffMassXiVsEffMassOmegaPlus);
 497 }
 498
 499 //-------
 500
 501 if(! f2dHistXiRadiusVsEffMassXiMinus) {
 502         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);
 503         fListHistCascade->Add(f2dHistXiRadiusVsEffMassXiMinus);
 504 }
 505
 506 if(! f2dHistXiRadiusVsEffMassXiPlus) {
 507         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);
 508         fListHistCascade->Add(f2dHistXiRadiusVsEffMassXiPlus);
 509 }
 510
 511 if(! f2dHistXiRadiusVsEffMassOmegaMinus) {
 512         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);
 513         fListHistCascade->Add(f2dHistXiRadiusVsEffMassOmegaMinus);
 514 }
 515
 516 if(! f2dHistXiRadiusVsEffMassOmegaPlus) {
 517         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);
 518         fListHistCascade->Add(f2dHistXiRadiusVsEffMassOmegaPlus);
 519 }
 520
 521
 522 //-------
 523
 524 if(! f2dHistXiPtVsEffMassXiMinus) {
 525         f2dHistXiPtVsEffMassXiMinus = new TH2F( "f2dHistXiPtVsEffMassXiMinus", "Pt_{cascade} Vs M_{#Xi^{-} candidates}; Pt_{cascade} (GeV/c); M( #Lambda , #pi^{-} ) (GeV/c^{2}) ", 100, 0., 10.0, 200, 1.2, 2.0);
 526         fListHistCascade->Add(f2dHistXiPtVsEffMassXiMinus);
 527 }
 528
 529 if(! f2dHistXiPtVsEffMassXiPlus) {
 530         f2dHistXiPtVsEffMassXiPlus = new TH2F( "f2dHistXiPtVsEffMassXiPlus", "Pt_{cascade} Vs M_{#Xi^{+} candidates}; Pt_{cascade} (GeV/c); M( #Lambda , #pi^{+} ) (GeV/c^{2}) ", 100, 0., 10.0, 200, 1.2, 2.0);
 531         fListHistCascade->Add(f2dHistXiPtVsEffMassXiPlus);
 532 }
 533
 534 if(! f2dHistXiPtVsEffMassOmegaMinus) {
 535         f2dHistXiPtVsEffMassOmegaMinus = new TH2F( "f2dHistXiPtVsEffMassOmegaMinus", "Pt_{cascade} Vs M_{#Omega^{-} candidates}; Pt_{cascade} (GeV/c); M( #Lambda , K^{-} ) (GeV/c^{2}) ", 100, 0., 10.0, 250, 1.5, 2.5);
 536         fListHistCascade->Add(f2dHistXiPtVsEffMassOmegaMinus);
 537 }
 538
 539 if(! f2dHistXiPtVsEffMassOmegaPlus) {
 540         f2dHistXiPtVsEffMassOmegaPlus = new TH2F( "f2dHistXiPtVsEffMassOmegaPlus", "Pt_{cascade} Vs M_{#Omega^{+} candidates}; Pt_{cascade} (GeV/c); M( #Lambda , K^{+} ) (GeV/c^{2}) ", 100, 0., 10.0, 250, 1.5, 2.5);
 541         fListHistCascade->Add(f2dHistXiPtVsEffMassOmegaPlus);
 542 }
 543
 544
 545 }// end UserCreateOutputObjects
 546
 547
 548
 549
 550
 551
 552 //________________________________________________________________________
 553 void AliAnalysisTaskCheckCascade::UserExec(Option_t *)
 554 {
 555   // Main loop
 556   // Called for each event
 557
 558         AliESDEvent *lESDevent = 0x0;
 559         AliAODEvent *lAODevent = 0x0;
 560         Int_t ncascades = -1;
 561
 562
 563   // Connect to the InputEvent
 564   // After these lines, we should have an ESD/AOD event + the number of cascades in it.
 565
 566   if(fAnalysisType == "ESD"){
 567         lESDevent = dynamic_cast<AliESDEvent*>( InputEvent() );
 568         if (!lESDevent) {
 569                 Printf("ERROR: lESDevent not available \n");
 570                 return;
 571         }
 572         ncascades = lESDevent->GetNumberOfCascades();
 573   }
 574
 575   if(fAnalysisType == "AOD"){
 576           lAODevent = dynamic_cast<AliAODEvent*>( InputEvent() );
 577         if (!lAODevent) {
 578                 Printf("ERROR: lAODevent not available \n");
 579                 return;
 580         }
 581         ncascades = lAODevent->GetNumberOfCascades();
 582   }
 583
 584   // ---------------------------------------------------------------
 585   // I - Initialisation of the local variables that will be needed
 586
 587
 588         // - 0th part of initialisation : around primary vertex ...
 589   Short_t  lStatusTrackingPrimVtx  = -2;
 590   Double_t lTrkgPrimaryVtxPos[3]   = {-100.0, -100.0, -100.0};
 591   Double_t lTrkgPrimaryVtxRadius3D = -500.0;
 592
 593   Double_t lBestPrimaryVtxPos[3]   = {-100.0, -100.0, -100.0};
 594   Double_t lBestPrimaryVtxRadius3D = -500.0;
 595
 596         // - 1st part of initialisation : variables needed to store AliESDCascade data members
 597   Double_t lEffMassXi      = 0. ;
 598   Double_t lChi2Xi         = 0. ;
 599   Double_t lDcaXiDaughters = 0. ;
 600   Double_t lXiCosineOfPointingAngle = 0. ;
 601   Double_t lPosXi[3] = { -1000.0, -1000.0, -1000.0 };
 602   Double_t lXiRadius = 0. ;
 603
 604
 605         // - 2nd part of initialisation : about V0 part in cascades
 606
 607   Double_t lInvMassLambdaAsCascDghter = 0.;
 608   Double_t lV0Chi2Xi         = 0. ;
 609   Double_t lDcaV0DaughtersXi = 0.;
 610
 611   Double_t lDcaBachToPrimVertexXi = 0., lDcaV0ToPrimVertexXi = 0.;
 612   Double_t lDcaPosToPrimVertexXi  = 0.;
 613   Double_t lDcaNegToPrimVertexXi  = 0.;
 614   Double_t lV0CosineOfPointingAngleXi = 0. ;
 615   Double_t lPosV0Xi[3] = { -1000. , -1000., -1000. }; // Position of VO coming from cascade
 616   Double_t lV0RadiusXi = -1000.0;
 617   Double_t lV0quality  = 0.;
 618
 619
 620         // - 3rd part of initialisation : Effective masses
 621
 622   Double_t lInvMassXiMinus    = 0.;
 623   Double_t lInvMassXiPlus     = 0.;
 624   Double_t lInvMassOmegaMinus = 0.;
 625   Double_t lInvMassOmegaPlus  = 0.;
 626
 627         // - 4th part of initialisation : PID treatment of the bachelor track
 628   Bool_t   lIsBachelorKaon    = kFALSE;
 629   Bool_t   lIsBachelorPion    = kFALSE;
 630
 631
 632         // - 5th part of initialisation : extra info for QA
 633
 634   Double_t lXiMomX       = 0., lXiMomY = 0., lXiMomZ = 0.;
 635   Double_t lXiTransvMom  = 0. ;
 636   Double_t lXiTotMom     = 0. ;
 637
 638   Double_t lBachMomX       = 0., lBachMomY  = 0., lBachMomZ   = 0.;
 639   Double_t lBachTransvMom  = 0.;
 640   Double_t lBachTotMom     = 0.;
 641
 642   Short_t  lChargeXi = -2;
 643   Double_t lV0toXiCosineOfPointingAngle = 0. ;
 644
 645   Double_t lRapXi   = -20.0, lRapOmega = -20.0,  lEta = -20.0, lTheta = 360., lPhi = 720. ;
 646   Double_t lAlphaXi = -200., lPtArmXi  = -200.0;
 647
 648
 649
 650   //-------------------------------------------------
 651   // O - Cascade vertexer (ESD)
 652
 653         //   if(fAnalysisType == "ESD" ){
 654         //      lESDevent->ResetCascades();
 655         //      AliCascadeVertexer CascVtxer;
 656         //      CascVtxer.V0sTracks2CascadeVertices(lESDevent);
 657         //     }
 658
 659
 660   // ---------------------------------------------------------------
 661   // I - General histos (filled for any event) - (ESD)
 662
 663   fHistTrackMultiplicity  ->Fill( (InputEvent())->GetNumberOfTracks() );
 664   fHistCascadeMultiplicity->Fill( ncascades );
 665
 666
 667
 668
 669   for (Int_t iXi = 0; iXi < ncascades; iXi++)
 670   {// This is the begining of the Cascade loop (ESD or AOD)
 671
 672
 673   if(fAnalysisType == "ESD"){
 674
 675   // -------------------------------------
 676   // II - Calcultaion Part dedicated to Xi vertices (ESD)
 677
 678         AliESDcascade *xi = lESDevent->GetCascade(iXi);
 679         if (!xi) continue;
 680
 681         // Just to know which file is currently open : locate the file containing Xi
 682         // cout << "Name of the file containing Xi candidate(s) :"
 683         //      << fInputHandler->GetTree()->GetCurrentFile()->GetName()
 684         //      << endl;
 685
 686
 687                 // - II.Step 1 : Characteristics of the event : prim. Vtx + magnetic field (ESD)
 688                 //-------------
 689
 690         // For new code (v4-16-Release-Rev06 or trunk)
 691
 692         const AliESDVertex *lPrimaryTrackingVtx = lESDevent->GetPrimaryVertexTracks();
 693                 // get the vtx stored in ESD found with tracks
 694
 695                 lPrimaryTrackingVtx->GetXYZ( lTrkgPrimaryVtxPos );
 696                 lTrkgPrimaryVtxRadius3D = TMath::Sqrt(  lTrkgPrimaryVtxPos[0] * lTrkgPrimaryVtxPos[0] +
 697                                                         lTrkgPrimaryVtxPos[1] * lTrkgPrimaryVtxPos[1] +
 698                                                         lTrkgPrimaryVtxPos[2] * lTrkgPrimaryVtxPos[2] );
 699
 700                 lStatusTrackingPrimVtx = lPrimaryTrackingVtx->GetStatus();
 701
 702
 703         const AliESDVertex *lPrimaryBestVtx = lESDevent->GetPrimaryVertex();
 704                 // get the best primary vertex available for the event
 705                 // As done in AliCascadeVertexer, we keep the one which is the best one available.
 706                 // between : Tracking vertex > SPD vertex > TPC vertex > default SPD vertex
 707                 // This one will be used for next calculations (DCA essentially)
 708
 709                 lPrimaryBestVtx->GetXYZ( lBestPrimaryVtxPos );
 710                 lBestPrimaryVtxRadius3D = TMath::Sqrt(  lBestPrimaryVtxPos[0] * lBestPrimaryVtxPos[0] +
 711                                                         lBestPrimaryVtxPos[1] * lBestPrimaryVtxPos[1] +
 712                                                         lBestPrimaryVtxPos[2] * lBestPrimaryVtxPos[2] );
 713
 714
 715         // For older evts
 716
 717         //      const AliESDVertex *lPrimaryTrackingVtx = lESDevent->GetPrimaryVertexTracks();
 718         //              get the vtx stored in ESD found with tracks
 719         //      Double_t lTrkgPrimaryVtxPos[3] = {-100.0, -100.0, -100.0};
 720         //              lPrimaryTrackingVtx->GetXYZ( lTrkgPrimaryVtxPos );
 721         //              Double_t lTrkgPrimaryVtxRadius3D = TMath::Sqrt( lTrkgPrimaryVtxPos[0]*lTrkgPrimaryVtxPos[0] +
 722         //                                              lTrkgPrimaryVtxPos[1] * lTrkgPrimaryVtxPos[1] +
 723         //                                              lTrkgPrimaryVtxPos[2] * lTrkgPrimaryVtxPos[2] );
 724         //
 725         //      const AliESDVertex *lPrimarySPDVtx = lESDevent->GetPrimaryVertexSPD();
 726         //              get the vtx found by exclusive use of SPD
 727         //      Double_t lSPDPrimaryVtxPos[3] = {-100.0, -100.0, -100.0};
 728         //              lPrimarySPDVtx->GetXYZ( lSPDPrimaryVtxPos );
 729         //
 730         //      // As done in AliCascadeVertexer, we keep, between both retrieved vertices,
 731         //      // the one which is the best one available.
 732         //      // This one will be used for next calculations (DCA essentially)
 733         //
 734         //              Double_t lBestPrimaryVtxPos[3] = {-100.0, -100.0, -100.0};
 735         //              if( lPrimaryTrackingVtx->GetStatus() ) { // if tracking vtx = ok
 736         //                      lBestPrimaryVtxPos[0] = lTrkgPrimaryVtxPos[0];
 737         //                      lBestPrimaryVtxPos[1] = lTrkgPrimaryVtxPos[1];
 738         //                      lBestPrimaryVtxPos[2] = lTrkgPrimaryVtxPos[2];
 739         //              }
 740         //              else{
 741         //                      lBestPrimaryVtxPos[0] = lSPDPrimaryVtxPos[0];
 742         //                      lBestPrimaryVtxPos[1] = lSPDPrimaryVtxPos[1];
 743         //                      lBestPrimaryVtxPos[2] = lSPDPrimaryVtxPos[2];
 744         //              }
 745         //
 746         //      Double_t lBestPrimaryVtxRadius3D = TMath::Sqrt( lBestPrimaryVtxPos[0]*lBestPrimaryVtxPos[0] +
 747         //                                              lBestPrimaryVtxPos[1] * lBestPrimaryVtxPos[1] +
 748         //                                              lBestPrimaryVtxPos[2] * lBestPrimaryVtxPos[2] );
 749
 750         // Back to normal
 751
 752         Double_t lMagneticField = lESDevent->GetMagneticField( );
 753
 754
 755
 756                 // - II.Step 2 : Assigning the necessary variables for specific AliESDcascade data members (ESD)
 757                 //-------------
 758         lV0quality = 0.;
 759         xi->ChangeMassHypothesis(lV0quality , 3312); // default working hypothesis : cascade = Xi- decay
 760
 761         lEffMassXi                      = xi->GetEffMassXi();
 762         lChi2Xi                         = xi->GetChi2Xi();
 763         lDcaXiDaughters                 = xi->GetDcaXiDaughters();
 764         lXiCosineOfPointingAngle        = xi->GetCascadeCosineOfPointingAngle( lBestPrimaryVtxPos[0],
 765                                                                                 lBestPrimaryVtxPos[1],
 766                                                                                 lBestPrimaryVtxPos[2] );
 767                 // Take care : the best available vertex should be used (like in AliCascadeVertexer)
 768
 769         xi->GetXYZcascade( lPosXi[0],  lPosXi[1], lPosXi[2] );
 770         lXiRadius                       = TMath::Sqrt( lPosXi[0]*lPosXi[0]  +  lPosXi[1]*lPosXi[1] );
 771
 772
 773
 774                 // - II.Step 3 : around the tracks : Bach + V0 (ESD)
 775                 // ~ Necessary variables for ESDcascade data members coming from the ESDv0 part (inheritance)
 776                 //-------------
 777
 778                 UInt_t lIdxPosXi        = (UInt_t) TMath::Abs( xi->GetPindex() );
 779                 UInt_t lIdxNegXi        = (UInt_t) TMath::Abs( xi->GetNindex() );
 780                 UInt_t lBachIdx         = (UInt_t) TMath::Abs( xi->GetBindex() );
 781                         // Care track label can be negative in MC production (linked with the track quality)
 782                         // However = normally, not the case for track index ...
 783
 784         AliESDtrack *pTrackXi           = lESDevent->GetTrack( lIdxPosXi );
 785         AliESDtrack *nTrackXi           = lESDevent->GetTrack( lIdxNegXi );
 786         AliESDtrack *bachTrackXi        = lESDevent->GetTrack( lBachIdx );
 787         if (!pTrackXi || !nTrackXi || !bachTrackXi ) {
 788                 Printf("ERROR: Could not retrieve one of the 3 daughter tracks of the cascade ...");
 789                 continue;
 790         }
 791
 792         lInvMassLambdaAsCascDghter      = xi->GetEffMass();
 793                 // This value shouldn't change, whatever the working hyp. is : Xi-, Xi+, Omega-, Omega+
 794         lDcaV0DaughtersXi               = xi->GetDcaV0Daughters();
 795         lV0Chi2Xi                       = xi->GetChi2V0();
 796
 797         lV0CosineOfPointingAngleXi      = xi->GetV0CosineOfPointingAngle( lBestPrimaryVtxPos[0],
 798                                                                           lBestPrimaryVtxPos[1],
 799                                                                           lBestPrimaryVtxPos[2] );
 800
 801         lDcaV0ToPrimVertexXi            = xi->GetD( lBestPrimaryVtxPos[0],
 802                                                     lBestPrimaryVtxPos[1],
 803                                                     lBestPrimaryVtxPos[2] );
 804
 805         lDcaBachToPrimVertexXi = TMath::Abs( bachTrackXi->GetD( lBestPrimaryVtxPos[0],
 806                                                                 lBestPrimaryVtxPos[1],
 807                                                                 lMagneticField  ) );
 808                                         // Note : AliExternalTrackParam::GetD returns an algebraic value ...
 809
 810                 xi->GetXYZ( lPosV0Xi[0],  lPosV0Xi[1], lPosV0Xi[2] );
 811         lV0RadiusXi             = TMath::Sqrt( lPosV0Xi[0]*lPosV0Xi[0]  +  lPosV0Xi[1]*lPosV0Xi[1] );
 812
 813         lDcaPosToPrimVertexXi   = TMath::Abs( pTrackXi  ->GetD( lBestPrimaryVtxPos[0],
 814                                                                 lBestPrimaryVtxPos[1],
 815                                                                 lMagneticField  )     );
 816
 817         lDcaNegToPrimVertexXi   = TMath::Abs( nTrackXi  ->GetD( lBestPrimaryVtxPos[0],
 818                                                                 lBestPrimaryVtxPos[1],
 819                                                                 lMagneticField  )     );
 820
 821
 822                 // - II.Step 4 : around effective masses (ESD)
 823                 // ~ change mass hypotheses to cover all the possibilities :  Xi-/+, Omega -/+
 824                 //-------------
 825
 826
 827         if( bachTrackXi->Charge() < 0 ) {
 828                 lV0quality = 0.;
 829                 xi->ChangeMassHypothesis(lV0quality , 3312);
 830                         // Calculate the effective mass of the Xi- candidate.
 831                         // pdg code 3312 = Xi-
 832                 lInvMassXiMinus = xi->GetEffMassXi();
 833
 834                 lV0quality = 0.;
 835                 xi->ChangeMassHypothesis(lV0quality , 3334);
 836                         // Calculate the effective mass of the Xi- candidate.
 837                         // pdg code 3334 = Omega-
 838                 lInvMassOmegaMinus = xi->GetEffMassXi();
 839
 840                 lV0quality = 0.;
 841                 xi->ChangeMassHypothesis(lV0quality , 3312);    // Back to default hyp.
 842         }// end if negative bachelor
 843
 844
 845         if( bachTrackXi->Charge() >  0 ){
 846                 lV0quality = 0.;
 847                 xi->ChangeMassHypothesis(lV0quality , -3312);
 848                         // Calculate the effective mass of the Xi+ candidate.
 849                         // pdg code -3312 = Xi+
 850                 lInvMassXiPlus = xi->GetEffMassXi();
 851
 852                 lV0quality = 0.;
 853                 xi->ChangeMassHypothesis(lV0quality , -3334);
 854                         // Calculate the effective mass of the Xi+ candidate.
 855                         // pdg code -3334  = Omega+
 856                 lInvMassOmegaPlus = xi->GetEffMassXi();
 857
 858                 lV0quality = 0.;
 859                 xi->ChangeMassHypothesis(lV0quality , -3312);   // Back to "default" hyp.
 860         }// end if positive bachelor
 861
 862
 863
 864                 // - II.Step 5 : PID on the bachelor
 865                 //-------------
 866
 867         // Reasonable guess for the priors for the cascade track sample
 868         Double_t lPriorsGuessXi[5]    = {0.0, 0.0, 2, 0, 1};
 869         Double_t lPriorsGuessOmega[5] = {0.0, 0.0, 1, 1, 1};
 870         AliPID pidXi;           pidXi.SetPriors(    lPriorsGuessXi    );
 871         AliPID pidOmega;        pidOmega.SetPriors( lPriorsGuessOmega );
 872
 873         if( bachTrackXi->IsOn(AliESDtrack::kESDpid) ){  // Combined PID exists
 874                 Double_t r[10]; bachTrackXi->GetESDpid(r);
 875                 pidXi.SetProbabilities(r);
 876                 pidOmega.SetProbabilities(r);
 877                 // Check if the bachelor track is a pion
 878                 Double_t ppion = pidXi.GetProbability(AliPID::kPion);
 879                 if (ppion > pidXi.GetProbability(AliPID::kElectron) &&
 880                     ppion > pidXi.GetProbability(AliPID::kMuon)     &&
 881                     ppion > pidXi.GetProbability(AliPID::kKaon)     &&
 882                     ppion > pidXi.GetProbability(AliPID::kProton)   )     lIsBachelorPion = kTRUE;
 883                 // Check if the bachelor track is a kaon
 884                 Double_t pkaon = pidOmega.GetProbability(AliPID::kKaon);
 885                 if (pkaon > pidOmega.GetProbability(AliPID::kElectron) &&
 886                     pkaon > pidOmega.GetProbability(AliPID::kMuon)     &&
 887                     pkaon > pidOmega.GetProbability(AliPID::kPion)     &&
 888                     pkaon > pidOmega.GetProbability(AliPID::kProton)   )  lIsBachelorKaon = kTRUE;
 889
 890         }// end if bachelor track with existing combined PID
 891
 892
 893
 894                 // - II.Step 6 : extra info for QA (ESD)
 895                 // miscellaneous pieces onf info that may help regarding data quality assessment.
 896                 //-------------
 897
 898         xi->GetPxPyPz( lXiMomX, lXiMomY, lXiMomZ );
 899                 lXiTransvMom    = TMath::Sqrt( lXiMomX*lXiMomX   + lXiMomY*lXiMomY );
 900                 lXiTotMom       = TMath::Sqrt( lXiMomX*lXiMomX   + lXiMomY*lXiMomY   + lXiMomZ*lXiMomZ );
 901
 902         xi->GetBPxPyPz(  lBachMomX,  lBachMomY,  lBachMomZ );
 903                 lBachTransvMom  = TMath::Sqrt( lBachMomX*lBachMomX   + lBachMomY*lBachMomY );
 904                 lBachTotMom     = TMath::Sqrt( lBachMomX*lBachMomX   + lBachMomY*lBachMomY  +  lBachMomZ*lBachMomZ  );
 905
 906         lChargeXi = xi->Charge();
 907
 908         lV0toXiCosineOfPointingAngle = xi->GetV0CosineOfPointingAngle( lPosXi[0], lPosXi[1], lPosXi[2] );
 909
 910         lRapXi    = xi->RapXi();
 911         lRapOmega = xi->RapOmega();
 912         lEta      = xi->Eta();
 913         lTheta    = xi->Theta() *180.0/TMath::Pi();
 914         lPhi      = xi->Phi()   *180.0/TMath::Pi();
 915         lAlphaXi  = xi->AlphaXi();
 916         lPtArmXi  = xi->PtArmXi();
 917
 918
 919   }// end of ESD treatment
 920
 921
 922   if(fAnalysisType == "AOD"){
 923
 924         // -------------------------------------
 925         // II - Calcultaion Part dedicated to Xi vertices (ESD)
 926
 927         const AliAODcascade *xi = lAODevent->GetCascade(iXi);
 928         if (!xi) continue;
 929
 930         // Just to know which file is currently open : locate the file containing Xi
 931         // cout << "Name of the file containing Xi candidate(s) :" <<  fesdH->GetTree()->GetCurrentFile()->GetName() << endl;
 932
 933
 934                 // - II.Step 1 : Characteristics of the event : prim. Vtx + magnetic field (AOD)
 935                 //-------------
 936
 937         lTrkgPrimaryVtxPos[0]   = -100.0;
 938         lTrkgPrimaryVtxPos[1]   = -100.0;
 939         lTrkgPrimaryVtxPos[2]   = -100.0;
 940         lTrkgPrimaryVtxRadius3D = -500. ;
 941         // We don't have the different prim. vertex at the AOD level -> nothing to do.
 942
 943         const AliAODVertex *lPrimaryBestVtx = lAODevent->GetPrimaryVertex();
 944         // get the best primary vertex available for the event
 945         // We may keep the one which is the best one available = GetVertex(0)
 946         // Pb with pile-up to expect
 947         // This one will be used for next calculations (DCA essentially)
 948
 949         lPrimaryBestVtx->GetXYZ( lBestPrimaryVtxPos );
 950         lBestPrimaryVtxRadius3D = TMath::Sqrt(  lBestPrimaryVtxPos[0] * lBestPrimaryVtxPos[0] +
 951                                                 lBestPrimaryVtxPos[1] * lBestPrimaryVtxPos[1] +
 952                                                 lBestPrimaryVtxPos[2] * lBestPrimaryVtxPos[2] );
 953
 954
 955                 // - II.Step 2 : Assigning the necessary variables for specific AliAODcascade data members (AOD)
 956                 //-------------
 957
 958         lEffMassXi                      = xi->MassXi(); // default working hypothesis : cascade = Xi- decay
 959         lChi2Xi                         = xi->Chi2Xi();
 960         lDcaXiDaughters                 = xi->DcaXiDaughters();
 961         lXiCosineOfPointingAngle        = xi->CosPointingAngleXi( lBestPrimaryVtxPos[0],
 962                                                                   lBestPrimaryVtxPos[1],
 963                                                                   lBestPrimaryVtxPos[2] );
 964                                         // Take care :
 965                                         // the best available vertex should be used (like in AliCascadeVertexer)
 966
 967                 lPosXi[0] = xi->DecayVertexXiX();
 968                 lPosXi[1] = xi->DecayVertexXiY();
 969                 lPosXi[2] = xi->DecayVertexXiZ();
 970         lXiRadius = TMath::Sqrt( lPosXi[0]*lPosXi[0]  +  lPosXi[1]*lPosXi[1] );
 971
 972
 973                 // - II.Step 3 : around the tracks : Bach + V0 (AOD)
 974                 // ~ Necessary variables for AODcascade data members coming from the AODv0 part (inheritance)
 975                 //-------------
 976
 977         lChargeXi                       = xi->ChargeXi();
 978
 979         if( lChargeXi < 0)
 980           lInvMassLambdaAsCascDghter    = xi->MassLambda();
 981         else
 982           lInvMassLambdaAsCascDghter    = xi->MassAntiLambda();
 983
 984         lDcaV0DaughtersXi               = xi->DcaV0Daughters();
 985         lV0Chi2Xi                       = xi->Chi2V0();
 986
 987         lV0CosineOfPointingAngleXi      = xi->CosPointingAngle( lBestPrimaryVtxPos );
 988         lDcaV0ToPrimVertexXi            = xi->DcaV0ToPrimVertex();
 989
 990         lDcaBachToPrimVertexXi          = xi->DcaBachToPrimVertex();
 991
 992
 993                 lPosV0Xi[0] = xi->DecayVertexV0X();
 994                 lPosV0Xi[1] = xi->DecayVertexV0Y();
 995                 lPosV0Xi[2] = xi->DecayVertexV0Z();
 996         lV0RadiusXi     = TMath::Sqrt( lPosV0Xi[0]*lPosV0Xi[0]  +  lPosV0Xi[1]*lPosV0Xi[1] );
 997
 998         lDcaPosToPrimVertexXi           = xi->DcaPosToPrimVertex();
 999         lDcaNegToPrimVertexXi           = xi->DcaNegToPrimVertex();
1000
1001
1002                 // - II.Step 4 : around effective masses (AOD)
1003                 // ~ change mass hypotheses to cover all the possibilities :  Xi-/+, Omega -/+
1004                 //-------------
1005
1006         if( lChargeXi < 0 )             lInvMassXiMinus         = xi->MassXi();
1007         if( lChargeXi > 0 )             lInvMassXiPlus          = xi->MassXi();
1008         if( lChargeXi < 0 )             lInvMassOmegaMinus      = xi->MassOmega();
1009         if( lChargeXi > 0 )             lInvMassOmegaPlus       = xi->MassOmega();
1010
1011
1012                 // - II.Step 5 : PID on the bachelor
1013                 //-------------
1014
1015
1016         // To be developed
1017
1018                 // - II.Step 6 : extra info for QA (AOD)
1019                 // miscellaneous pieces onf info that may help regarding data quality assessment.
1020                 //-------------
1021
1022                 lXiMomX = xi->MomXiX();
1023                 lXiMomY = xi->MomXiY();
1024                 lXiMomZ = xi->MomXiZ();
1025         lXiTransvMom    = TMath::Sqrt( lXiMomX*lXiMomX   + lXiMomY*lXiMomY );
1026         lXiTotMom       = TMath::Sqrt( lXiMomX*lXiMomX   + lXiMomY*lXiMomY   + lXiMomZ*lXiMomZ );
1027
1028                 lBachMomX = xi->MomBachX();
1029                 lBachMomY = xi->MomBachY();
1030                 lBachMomZ = xi->MomBachZ();
1031         lBachTransvMom  = TMath::Sqrt( lBachMomX*lBachMomX   + lBachMomY*lBachMomY );
1032         lBachTotMom     = TMath::Sqrt( lBachMomX*lBachMomX   + lBachMomY*lBachMomY  +  lBachMomZ*lBachMomZ  );
1033
1034
1035         lV0toXiCosineOfPointingAngle = xi->CosPointingAngle( xi->GetDecayVertexXi() );
1036
1037         lRapXi    = xi->RapXi();
1038         lRapOmega = xi->RapOmega();
1039         lEta      = xi->Eta();                          // Will not work ! need a method Pz(), Py() Px()
1040         lTheta    = xi->Theta() *180.0/TMath::Pi();     // in AODcascade.
1041         lPhi      = xi->Phi()   *180.0/TMath::Pi();     // Here, we will get eta, theta, phi for the V0 ...
1042         lAlphaXi  = xi->AlphaXi();
1043         lPtArmXi  = xi->PtArmXi();
1044
1045
1046   }// end of AOD treatment
1047
1048
1049   // -------------------------------------
1050   // III - Filling the TH1Fs
1051
1052
1053         // - III.Step 1
1054
1055         fHistVtxStatus                  ->Fill( lStatusTrackingPrimVtx   );  // 1 if tracking vtx = ok
1056
1057         if( lStatusTrackingPrimVtx ){
1058                 fHistPosTrkgPrimaryVtxX ->Fill( lTrkgPrimaryVtxPos[0]    );
1059                 fHistPosTrkgPrimaryVtxY ->Fill( lTrkgPrimaryVtxPos[1]    );
1060                 fHistPosTrkgPrimaryVtxZ ->Fill( lTrkgPrimaryVtxPos[2]    );
1061                 fHistTrkgPrimaryVtxRadius->Fill( lTrkgPrimaryVtxRadius3D );
1062         }
1063
1064         fHistPosBestPrimaryVtxX         ->Fill( lBestPrimaryVtxPos[0]    );
1065         fHistPosBestPrimaryVtxY         ->Fill( lBestPrimaryVtxPos[1]    );
1066         fHistPosBestPrimaryVtxZ         ->Fill( lBestPrimaryVtxPos[2]    );
1067         fHistBestPrimaryVtxRadius       ->Fill( lBestPrimaryVtxRadius3D  );
1068
1069         f2dHistTrkgPrimVtxVsBestPrimVtx->Fill( lTrkgPrimaryVtxRadius3D, lBestPrimaryVtxRadius3D );
1070
1071
1072         // - III.Step 2
1073         fHistEffMassXi                  ->Fill( lEffMassXi               );
1074         fHistChi2Xi                     ->Fill( lChi2Xi                  );     // Flag CascadeVtxer: Cut Variable a
1075         fHistDcaXiDaughters             ->Fill( lDcaXiDaughters          );     // Flag CascadeVtxer: Cut Variable e
1076         fHistDcaBachToPrimVertex        ->Fill( lDcaBachToPrimVertexXi   );     // Flag CascadeVtxer: Cut Variable d
1077         fHistXiCosineOfPointingAngle    ->Fill( lXiCosineOfPointingAngle );     // Flag CascadeVtxer: Cut Variable f
1078         fHistXiRadius                   ->Fill( lXiRadius                );     // Flag CascadeVtxer: Cut Variable g+h
1079
1080
1081         // - III.Step 3
1082         fHistMassLambdaAsCascDghter     ->Fill( lInvMassLambdaAsCascDghter );   // Flag CascadeVtxer: Cut Variable c
1083         fHistV0Chi2Xi                   ->Fill( lV0Chi2Xi                  );
1084         fHistDcaV0DaughtersXi           ->Fill( lDcaV0DaughtersXi          );
1085         fHistV0CosineOfPointingAngleXi  ->Fill( lV0CosineOfPointingAngleXi );
1086         fHistV0RadiusXi                 ->Fill( lV0RadiusXi                );
1087
1088         fHistDcaV0ToPrimVertexXi        ->Fill( lDcaV0ToPrimVertexXi       );   // Flag CascadeVtxer: Cut Variable b
1089         fHistDcaPosToPrimVertexXi       ->Fill( lDcaPosToPrimVertexXi      );
1090         fHistDcaNegToPrimVertexXi       ->Fill( lDcaNegToPrimVertexXi      );
1091
1092
1093         // - III.Step 4+5
1094         if( lChargeXi < 0 ){
1095                                         fHistMassXiMinus               ->Fill( lInvMassXiMinus    );
1096                                         fHistMassOmegaMinus            ->Fill( lInvMassOmegaMinus );
1097                 if(lIsBachelorPion)     fHistMassWithCombPIDXiMinus    ->Fill( lInvMassXiMinus    );
1098                 if(lIsBachelorKaon)     fHistMassWithCombPIDOmegaMinus ->Fill( lInvMassOmegaMinus );
1099         }
1100
1101         if( lChargeXi > 0 ){
1102                                         fHistMassXiPlus                ->Fill( lInvMassXiPlus     );
1103                                         fHistMassOmegaPlus             ->Fill( lInvMassOmegaPlus  );
1104                 if(lIsBachelorPion)     fHistMassWithCombPIDXiPlus     ->Fill( lInvMassXiPlus     );
1105                 if(lIsBachelorKaon)     fHistMassWithCombPIDOmegaPlus  ->Fill( lInvMassOmegaPlus  );
1106         }
1107
1108
1109 //      if( bachTrackXi->Charge() < 0 ) fHistMassXiMinus        ->Fill( lInvMassXiMinus    );
1110 //      if( bachTrackXi->Charge() > 0 ) fHistMassXiPlus         ->Fill( lInvMassXiPlus     );
1111 //      if( bachTrackXi->Charge() < 0 ) fHistMassOmegaMinus     ->Fill( lInvMassOmegaMinus );
1112 //      if( bachTrackXi->Charge() > 0 ) fHistMassOmegaPlus      ->Fill( lInvMassOmegaPlus  );
1113
1114
1115         // - III.Step 5
1116         fHistXiTransvMom        ->Fill( lXiTransvMom   );
1117         fHistXiTotMom           ->Fill( lXiTotMom      );
1118
1119         fHistBachTransvMom      ->Fill( lBachTransvMom );
1120         fHistBachTotMom         ->Fill( lBachTotMom    );
1121
1122         fHistChargeXi           ->Fill( lChargeXi      );
1123         fHistV0toXiCosineOfPointingAngle->Fill( lV0toXiCosineOfPointingAngle );
1124
1125         fHistRapXi              ->Fill( lRapXi               );
1126         fHistRapOmega           ->Fill( lRapOmega            );
1127         fHistEta                ->Fill( lEta                 );
1128         fHistTheta              ->Fill( lTheta               );
1129         fHistPhi                ->Fill( lPhi                 );
1130
1131         f2dHistArmenteros       ->Fill( lAlphaXi, lPtArmXi   );
1132
1133         if( lChargeXi < 0 ) {
1134                 f2dHistEffMassLambdaVsEffMassXiMinus->Fill( lInvMassLambdaAsCascDghter, lInvMassXiMinus );
1135                 f2dHistEffMassXiVsEffMassOmegaMinus ->Fill( lInvMassXiMinus, lInvMassOmegaMinus );
1136                 f2dHistXiRadiusVsEffMassXiMinus     ->Fill( lXiRadius, lInvMassXiMinus );
1137                 f2dHistXiRadiusVsEffMassOmegaMinus  ->Fill( lXiRadius, lInvMassOmegaMinus );
1138                 f2dHistXiPtVsEffMassXiMinus         ->Fill( lXiTransvMom, lInvMassXiMinus );
1139                 f2dHistXiPtVsEffMassOmegaMinus      ->Fill( lXiTransvMom, lInvMassOmegaMinus );
1140         }
1141         else{
1142                 f2dHistEffMassLambdaVsEffMassXiPlus ->Fill( lInvMassLambdaAsCascDghter, lInvMassXiPlus );
1143                 f2dHistEffMassXiVsEffMassOmegaPlus  ->Fill( lInvMassXiPlus, lInvMassOmegaPlus );
1144                 f2dHistXiRadiusVsEffMassXiPlus      ->Fill( lXiRadius, lInvMassXiPlus);
1145                 f2dHistXiRadiusVsEffMassOmegaPlus   ->Fill( lXiRadius, lInvMassOmegaPlus );
1146                 f2dHistXiPtVsEffMassXiPlus          ->Fill( lXiTransvMom, lInvMassXiPlus );
1147                 f2dHistXiPtVsEffMassOmegaPlus       ->Fill( lXiTransvMom, lInvMassOmegaPlus );
1148         }
1149
1150
1151     }// end of the Cascade loop (ESD or AOD)
1152
1153
1154   // Post output data.
1155  PostData(1, fListHistCascade);
1156 }
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167 //________________________________________________________________________
1168 void AliAnalysisTaskCheckCascade::Terminate(Option_t *)
1169 {
1170   // Draw result to the screen
1171   // Called once at the end of the query
1172
1173   TList *cRetrievedList = 0x0;
1174          cRetrievedList = (TList*)GetOutputData(1);
1175         if(!cRetrievedList){
1176                 Printf("ERROR - AliAnalysisTaskCheckCascade: ouput data container list not available\n");
1177                 return;
1178         }
1179
1180   fHistTrackMultiplicity = dynamic_cast<TH1F*> (   cRetrievedList->FindObject("fHistTrackMultiplicity") );
1181         if (!fHistTrackMultiplicity) {
1182                 Printf("ERROR - AliAnalysisTaskCheckCascade: fHistTrackMultiplicity not available\n");
1183                 return;
1184         }
1185
1186   fHistCascadeMultiplicity = dynamic_cast<TH1F*> ( cRetrievedList->FindObject("fHistCascadeMultiplicity"));
1187         if (!fHistCascadeMultiplicity) {
1188                 Printf("ERROR - AliAnalysisTaskCheckCascade: fHistCascadeMultiplicity not available\n");
1189                 return;
1190         }
1191
1192   TCanvas *canCheckCascade = new TCanvas("AliAnalysisTaskCheckCascade","Multiplicity",10,10,510,510);
1193   canCheckCascade->cd(1)->SetLogy();
1194
1195   fHistTrackMultiplicity->SetMarkerStyle(22);
1196   fHistTrackMultiplicity->DrawCopy("E");
1197
1198   fHistCascadeMultiplicity->SetMarkerStyle(26);
1199   fHistCascadeMultiplicity->DrawCopy("ESAME");
1200
1201 }