PWGLF/RESONANCES/extra/AliXiStar.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 ////////////////////////////////////////////////////////////////////////////////
  17 //
  18 //  This class is used to reconstruct the neutral Xi(1530) resonance.
  19 //  This class essentially combines charged Xi candidates from the Xi Vertexer
  20 //  with primary charged pions.
  21 //
  22 //  authors: Dhevan Gangadharan (dhevan.raja.gangadharan@cern.ch)
  23 //
  24 ////////////////////////////////////////////////////////////////////////////////
  25
  26
  27
  28 #include <iostream>
  29 #include <math.h>
  30 #include "TChain.h"
  31 #include "TFile.h"
  32 #include "TKey.h"
  33 #include "TObject.h"
  34 #include "TObjString.h"
  35 #include "TList.h"
  36 #include "TTree.h"
  37 #include "TH1F.h"
  38 #include "TH1D.h"
  39 #include "TH2D.h"
  40 #include "TH3D.h"
  41 #include "TProfile.h"
  42 #include "TCanvas.h"
  43
  44 #include "AliAnalysisTask.h"
  45 #include "AliAnalysisManager.h"
  46
  47
  48 #include "AliESDEvent.h"
  49 #include "AliESDInputHandler.h"
  50 #include "AliESDtrackCuts.h"
  51 #include "AliMCEventHandler.h"
  52 #include "AliMCEvent.h"
  53 #include "AliStack.h"
  54
  55 #include "AliAODEvent.h"
  56 #include "AliAODInputHandler.h"
  57 #include "AliAODMCParticle.h"
  58 #include "AliAODcascade.h"
  59 #include "AliESDcascade.h"
  60 #include "AliV0vertexer.h"
  61 #include "AliCascadeVertexer.h"
  62
  63 #include "AliXiStar.h"
  64
  65 #define PI 3.1415927
  66
  67
  68 // Author: Dhevan Gangadharan
  69
  70 ClassImp(AliXiStar)
  71
  72 //________________________________________________________________________
  73 AliXiStar::AliXiStar():
  74 AliAnalysisTaskSE(),
  75   fname(0),
  76   fAOD(0x0),
  77   fESD(0x0),
  78   fOutputList(0x0),
  79   fTrackCut(0x0),
  80   fPIDResponse(0x0),
  81   fEC(0x0),
  82   fEvt(0x0),
  83   fTempStruct(0x0),
  84   fZvertexBins(0),
  85   fEventsToMix(0),
  86   fMultBins(0),
  87   fMCcase(0),
  88   fAODcase(0),
  89   fEventCounter(0),
  90   fMaxDecayLength(0),
  91   fMassWindow(0),
  92   fTrueMassPr(0),
  93   fTrueMassPi(0),
  94   fTrueMassK(0),
  95   fTrueMassLam(0),
  96   fTrueMassXi(0),
  97   fESDTrack4(0x0),
  98   fXiTrack(0x0),
  99   fCutList(0)
 100 {
 101   // Default Constructor
 102   for (Int_t i=0; i<21; i++){
 103     fCovMatrix[i]=-99999.;
 104     if (i<12) fMultLimits[i] = 0;
 105   }
 106   for (Int_t i=0; i<kNCuts; i++){
 107     fDecayParameters[i]=0;
 108     for (Int_t j=0; j<kNCutVariations; j++){
 109       fCutValues[j][i]=0;
 110     }
 111   }
 112   //
 113   for (Int_t cv=0; cv<kNCutVariations; cv++){
 114     CutVar[cv].fXi=0x0;
 115     CutVar[cv].fXibar=0x0;
 116     CutVar[cv].fXiMinusPiPlus=0x0;
 117     CutVar[cv].fXiMinusPiMinus=0x0;
 118     CutVar[cv].fXiPlusPiPlus=0x0;
 119     CutVar[cv].fXiPlusPiMinus=0x0;
 120     //
 121     CutVar[cv].fXiMinusPiPlusbkg=0x0;
 122     CutVar[cv].fXiMinusPiMinusbkg=0x0;
 123     CutVar[cv].fXiPlusPiPlusbkg=0x0;
 124     CutVar[cv].fXiPlusPiMinusbkg=0x0;
 125     //
 126     CutVar[cv].fMCrecXi=0x0;
 127     CutVar[cv].fMCrecXibar=0x0;
 128     CutVar[cv].fMCrecXiMinusPiPlus=0x0;
 129     CutVar[cv].fMCrecXiPlusPiMinus=0x0;
 130   }
 131
 132 }
 133 //________________________________________________________________________
 134 AliXiStar::AliXiStar(const char *name, Bool_t AODdecision, Bool_t MCdecision, Int_t CutListOption)
 135   : AliAnalysisTaskSE(name),
 136     fname(name),
 137     fAOD(0x0),
 138     fESD(0x0),
 139     fOutputList(0x0),
 140     fTrackCut(0x0),
 141     fPIDResponse(0x0),
 142     fEC(0x0),
 143     fEvt(0x0),
 144     fTempStruct(0x0),
 145     fZvertexBins(0),
 146     fEventsToMix(0),
 147     fMultBins(0),
 148     fMCcase(MCdecision),
 149     fAODcase(AODdecision),
 150     fEventCounter(0),
 151     fMaxDecayLength(0),
 152     fMassWindow(0),
 153     fTrueMassPr(0),
 154     fTrueMassPi(0),
 155     fTrueMassK(0),
 156     fTrueMassLam(0),
 157     fTrueMassXi(0),
 158     fESDTrack4(0x0),
 159     fXiTrack(0x0),
 160     fCutList(CutListOption)
 161
 162 {
 163   // Main Constructor
 164   for (Int_t i=0; i<21; i++){
 165     fCovMatrix[i]=-99999.;
 166     if (i<12) fMultLimits[i] = 0;
 167   }
 168   for (Int_t i=0; i<kNCuts; i++){
 169     fDecayParameters[i]=0;
 170     for (Int_t j=0; j<kNCutVariations; j++){
 171       fCutValues[j][i]=0;
 172     }
 173   }
 174   //
 175   for (Int_t cv=0; cv<kNCutVariations; cv++){
 176     CutVar[cv].fXi=0x0;
 177     CutVar[cv].fXibar=0x0;
 178     CutVar[cv].fXiMinusPiPlus=0x0;
 179     CutVar[cv].fXiMinusPiMinus=0x0;
 180     CutVar[cv].fXiPlusPiPlus=0x0;
 181     CutVar[cv].fXiPlusPiMinus=0x0;
 182     //
 183     CutVar[cv].fXiMinusPiPlusbkg=0x0;
 184     CutVar[cv].fXiMinusPiMinusbkg=0x0;
 185     CutVar[cv].fXiPlusPiPlusbkg=0x0;
 186     CutVar[cv].fXiPlusPiMinusbkg=0x0;
 187     //
 188     CutVar[cv].fMCrecXi=0x0;
 189     CutVar[cv].fMCrecXibar=0x0;
 190     CutVar[cv].fMCrecXiMinusPiPlus=0x0;
 191     CutVar[cv].fMCrecXiPlusPiMinus=0x0;
 192   }
 193
 194
 195   // Define output slots here
 196   // Output slot #1
 197   DefineOutput(1, TList::Class());
 198
 199 }
 200 //________________________________________________________________________
 201 AliXiStar::AliXiStar(const AliXiStar &obj)
 202   : AliAnalysisTaskSE(obj.fname),
 203     fname(obj.fname),
 204     fAOD(obj.fAOD),
 205     fESD(obj.fESD),
 206     fOutputList(obj.fOutputList),
 207     fTrackCut(obj.fTrackCut),
 208     fPIDResponse(obj.fPIDResponse),
 209     fEC(obj.fEC),
 210     fEvt(obj.fEvt),
 211     fTempStruct(obj.fTempStruct),
 212     fZvertexBins(obj.fZvertexBins),
 213     fEventsToMix(obj.fEventsToMix),
 214     fMultBins(obj.fMultBins),
 215     fMultLimits(),
 216     fMCcase(obj.fMCcase),
 217     fAODcase(obj.fAODcase),
 218     fEventCounter(obj.fEventCounter),
 219     fMaxDecayLength(obj.fMaxDecayLength),
 220     fMassWindow(obj.fMassWindow),
 221     fTrueMassPr(obj.fTrueMassPr),
 222     fTrueMassPi(obj.fTrueMassPi),
 223     fTrueMassK(obj.fTrueMassK),
 224     fTrueMassLam(obj.fTrueMassLam),
 225     fTrueMassXi(obj.fTrueMassXi),
 226     fESDTrack4(obj.fESDTrack4),
 227     fXiTrack(obj.fXiTrack),
 228     fCutList(obj.fCutList)
 229
 230 {
 231   // Copy constructor
 232   for (Int_t i=0; i<21; i++){
 233     fCovMatrix[i]=obj.fCovMatrix[i];
 234     if (i<12) fMultLimits[i]=obj.fMultLimits[i];
 235   }
 236   for (Int_t i=0; i<kNCuts; i++){
 237     fDecayParameters[i]=obj.fDecayParameters[i];
 238     for (Int_t j=0; j<kNCutVariations; j++){
 239       fCutValues[j][i]=obj.fCutValues[j][i];
 240     }
 241   }
 242
 243 }
 244 //________________________________________________________________________
 245 AliXiStar &AliXiStar::operator=(const AliXiStar &obj)
 246 {
 247   // Assignment operator
 248   if (this == &obj)
 249     return *this;
 250
 251   fname = obj.fname;
 252   fAOD = obj.fAOD;
 253   fESD = obj.fESD;
 254   fOutputList = obj.fOutputList;
 255   fTrackCut = obj.fTrackCut;
 256   fPIDResponse = obj.fPIDResponse;
 257   fEC = obj.fEC;
 258   fEvt = obj.fEvt;
 259   fTempStruct = obj.fTempStruct;
 260   fZvertexBins = obj.fZvertexBins;
 261   fEventsToMix = obj.fEventsToMix;
 262   fMultBins = obj.fMultBins;
 263   for (Int_t i=0; i<12; i++){
 264     fMultLimits[i]=obj.fMultLimits[i];
 265   }
 266   fMCcase = obj.fMCcase;
 267   fAODcase = obj.fAODcase;
 268   fEventCounter = obj.fEventCounter;
 269   fMaxDecayLength = obj.fMaxDecayLength;
 270   fMassWindow = obj.fMassWindow;
 271   for (Int_t i=0; i<21; i++){
 272     fCovMatrix[i]=obj.fCovMatrix[i];
 273   }
 274   fTrueMassPr = obj.fTrueMassPr;
 275   fTrueMassPi = obj.fTrueMassPi;
 276   fTrueMassK = obj.fTrueMassK;
 277   fTrueMassLam = obj.fTrueMassLam;
 278   fTrueMassXi = obj.fTrueMassXi;
 279   fESDTrack4 = obj.fESDTrack4;
 280   fXiTrack = obj.fXiTrack;
 281   fCutList = obj.fCutList;
 282
 283   for (Int_t i=0; i<kNCuts; i++){
 284     fDecayParameters[i]=obj.fDecayParameters[i];
 285     for (Int_t j=0; j<kNCutVariations; j++){
 286       fCutValues[j][i]=obj.fCutValues[j][i];
 287     }
 288   }
 289
 290
 291   return (*this);
 292 }
 293 //________________________________________________________________________
 294 AliXiStar::~AliXiStar()
 295 {
 296   // Destructor
 297   if(fAOD) delete fAOD;
 298   if(fESD) delete fESD;
 299   if(fOutputList) delete fOutputList;
 300   if(fTrackCut) delete fTrackCut;
 301   if(fPIDResponse) delete fPIDResponse;
 302   if(fEC) delete fEC;
 303   if(fEvt) delete fEvt;
 304   if(fTempStruct) delete fTempStruct;
 305   if(fESDTrack4) delete fESDTrack4;
 306   if(fXiTrack) delete fXiTrack;
 307
 308   for (Int_t cv=0; cv<kNCutVariations; cv++){
 309     if(CutVar[cv].fXi) delete CutVar[cv].fXi;
 310     if(CutVar[cv].fXibar) delete CutVar[cv].fXibar;
 311     if(CutVar[cv].fXiMinusPiPlus) delete CutVar[cv].fXiMinusPiPlus;
 312     if(CutVar[cv].fXiMinusPiMinus) delete CutVar[cv].fXiMinusPiMinus;
 313     if(CutVar[cv].fXiPlusPiPlus) delete CutVar[cv].fXiPlusPiPlus;
 314     if(CutVar[cv].fXiPlusPiMinus) delete CutVar[cv].fXiPlusPiMinus;
 315     //
 316     if(CutVar[cv].fXiMinusPiPlusbkg) delete CutVar[cv].fXiMinusPiPlusbkg;
 317     if(CutVar[cv].fXiMinusPiMinusbkg) delete CutVar[cv].fXiMinusPiMinusbkg;
 318     if(CutVar[cv].fXiPlusPiPlusbkg) delete CutVar[cv].fXiPlusPiPlusbkg;
 319     if(CutVar[cv].fXiPlusPiMinusbkg) delete CutVar[cv].fXiPlusPiMinusbkg;
 320     //
 321     if(CutVar[cv].fMCrecXi) delete CutVar[cv].fMCrecXi;
 322     if(CutVar[cv].fMCrecXibar) delete CutVar[cv].fMCrecXibar;
 323     if(CutVar[cv].fMCrecXiMinusPiPlus) delete CutVar[cv].fMCrecXiMinusPiPlus;
 324     if(CutVar[cv].fMCrecXiPlusPiMinus) delete CutVar[cv].fMCrecXiPlusPiMinus;
 325   }
 326
 327 }
 328 //________________________________________________________________________
 329 void AliXiStar::XiStarInit()
 330 {
 331   //
 332   //Inits cuts and analysis settings
 333   //
 334
 335   fEventCounter=0;// event counter initialization
 336   cout<<"AliXiStar XiStarInit() call"<<endl;
 337
 338
 339   ///////////////////////////////////////////////
 340   // Track Cuts for ESD analysis
 341   fTrackCut = new AliESDtrackCuts();
 342   fTrackCut->SetPtRange(.15,1000);
 343   fTrackCut->SetAcceptKinkDaughters(kFALSE);
 344   //fTrackCut->SetMinNClustersTPC(70);
 345   fTrackCut->SetRequireTPCRefit(kTRUE);
 346   ////////////////////////////////////////////////
 347
 348   fZvertexBins = 20;
 349   fMultBins = 11;// This must also be set in AliXiStar.h
 350   if(fMCcase) fEventsToMix = 0;
 351   else fEventsToMix = 40;
 352
 353   // multiplicity edges for event mixing bins
 354   fMultLimits[0]=0, fMultLimits[1]=5, fMultLimits[2]=10, fMultLimits[3]=15, fMultLimits[4]=20, fMultLimits[5]=25;
 355   fMultLimits[6]=30, fMultLimits[7]=35, fMultLimits[8]=40, fMultLimits[9]=45, fMultLimits[10]=50, fMultLimits[11]=150;
 356
 357
 358   fEC = new AliXiStarEventCollection **[fZvertexBins];
 359   for(unsigned short i=0; i<fZvertexBins; i++){
 360
 361     fEC[i] = new AliXiStarEventCollection *[fMultBins];
 362
 363     for(unsigned short j=0; j<fMultBins; j++){
 364
 365       fEC[i][j] = new AliXiStarEventCollection(fEventsToMix+1);
 366     }
 367   }
 368
 369
 370   fTempStruct = new AliXiStarTrackStruct[kNbinsM];
 371
 372   fESDTrack4 = new AliESDtrack();
 373   fXiTrack = new AliESDtrack();
 374
 375
 376   fMaxDecayLength = 100.;
 377   fMassWindow = 0.006;
 378
 379   /////////////////////////////////////////////////////////////////////////////////////
 380   ///////////////////////
 381   // DecayParameters Key (number represents array index)
 382   // NclustersTPC: 0=proton, 1=pion first, 2=pion second, 3=pion third
 383   // DCAVtx: 4=proton, 5=pion first, 6=pion second, 7=lambda, 8=pion third
 384   // 9 = DCA proton-pion
 385   // 10 = DCA Lambda-pion
 386   // 11 = Rxy Lambda
 387   // 12 = Rxy Xi
 388   // 13 = Cos PA Lambda
 389   // 14 = Cos PA Xi
 390
 391   // Set Standard Reconstruction cut values
 392   fCutValues[0][0] = 70; fCutValues[0][1] = 70; fCutValues[0][2] = 70; fCutValues[0][3] = 70;
 393   fCutValues[0][4] = 0.04; fCutValues[0][5] = 0.04; fCutValues[0][6] = 0.05; fCutValues[0][7] = 0.07; fCutValues[0][8] = 2.0;
 394   fCutValues[0][9] = 1.6;
 395   fCutValues[0][10] = 1.6;
 396   fCutValues[0][11] = 1.4;
 397   fCutValues[0][12] = 0.8;
 398   fCutValues[0][13] = 0.97;
 399   fCutValues[0][14] = 0.97;
 400   for(int cv=1; cv<kNCutVariations; cv++){
 401     for(int ct=0; ct<kNCuts; ct++){
 402       fCutValues[cv][ct] = fCutValues[0][ct];
 403     }
 404   }
 405   // Set Systematic Variations
 406   fCutValues[1][0] = 80; fCutValues[1][1] = 80; fCutValues[1][2] = 80; fCutValues[1][3] = 80;// 80
 407   fCutValues[2][4] = 0.104;// 0.104
 408   fCutValues[3][5] = 0.104;// 0.104
 409   fCutValues[4][6] = 0.08;// 0.08
 410   fCutValues[5][7] = 0.1;// 0.1
 411   fCutValues[6][8] = 1.0;// 1.0
 412   fCutValues[7][9] = 0.94;// 0.94
 413   fCutValues[8][10] = 1.41;// 1.41
 414   fCutValues[9][11] = 4.39;// 4.39
 415   fCutValues[10][12] = 0.95;// 0.95
 416   fCutValues[11][13] = 0.99;// 0.99
 417   fCutValues[12][14] = 0.985;// 0.085
 418
 419
 420
 421
 422
 423   // PDG mass values
 424   fTrueMassPr=.93827, fTrueMassPi=.13957, fTrueMassK=.493677, fTrueMassLam=1.11568, fTrueMassXi=1.32171;
 425
 426   // The following CovMatrix is set so that PropogateToDCA() ignores track errors. Only used to propagate Xi to third pion for XiStar reconstruction
 427   for(Int_t i=0; i<21; i++) fCovMatrix[i]=0;
 428   fCovMatrix[0]=1, fCovMatrix[2]=1, fCovMatrix[5]=1, fCovMatrix[9]=1, fCovMatrix[14]=1, fCovMatrix[20]=1;
 429
 430
 431 }
 432 //________________________________________________________________________
 433 void AliXiStar::UserCreateOutputObjects()
 434 {
 435
 436   XiStarInit();// Initialize settings
 437
 438   // Create histograms
 439   fOutputList = new TList();
 440   fOutputList->SetOwner();
 441
 442   TH3F *fVertexDist1 = new TH3F("fVertexDist1","Vertex Distribution",20,-1,1, 20,-1,1, 600,-30,30);
 443   fVertexDist1->GetXaxis()->SetTitle("X Vertex (cm)");
 444   fVertexDist1->GetYaxis()->SetTitle("Y Vertex (cm)");
 445   fVertexDist1->GetZaxis()->SetTitle("Z Vertex (cm)");
 446   fOutputList->Add(fVertexDist1);
 447
 448   TH3F *fVertexDist3 = new TH3F("fVertexDist3","Vertex Distribution",20,-1,1, 20,-1,1, 600,-30,30);
 449   fVertexDist3->GetXaxis()->SetTitle("X Vertex (cm)");
 450   fVertexDist3->GetYaxis()->SetTitle("Y Vertex (cm)");
 451   fVertexDist3->GetZaxis()->SetTitle("Z Vertex (cm)");
 452   fOutputList->Add(fVertexDist3);
 453
 454   TH2F *fDCADist = new TH2F("fDCADist","DCA distribution",kNbinsM,-.5,kNbinsM-.5, 100,0,10);
 455   fOutputList->Add(fDCADist);
 456
 457
 458   TH3F *fMultDist3d = new TH3F("fMultDist3d","Multiplicity Distribution",kNbinsM,-.5,kNbinsM-.5, kNbinsM,-.5,kNbinsM-.5, kNbinsM,-.5,kNbinsM-.5);
 459   fMultDist3d->GetXaxis()->SetTitle("Multiplicity");
 460   fMultDist3d->GetYaxis()->SetTitle("Positive Multiplicity");
 461   fMultDist3d->GetZaxis()->SetTitle("Negative Multiplicity");
 462   fMultDist3d->SetMarkerStyle(kFullCircle);
 463   fOutputList->Add(fMultDist3d);
 464
 465
 466   TH1F *fMultDist1 = new TH1F("fMultDist1","Multiplicity Distribution",kNbinsM,-.5,kNbinsM-.5);
 467   fMultDist1->GetXaxis()->SetTitle("Multiplicity");
 468   fOutputList->Add(fMultDist1);
 469
 470   TH1F *fMultDist2 = new TH1F("fMultDist2","Multiplicity Distribution",kNbinsM,-.5,kNbinsM-.5);
 471   fMultDist2->GetXaxis()->SetTitle("Multiplicity");
 472   fOutputList->Add(fMultDist2);
 473
 474   TH1F *fMultDist3 = new TH1F("fMultDist3","Multiplicity Distribution",kNbinsM,-.5,kNbinsM-.5);
 475   fMultDist3->GetXaxis()->SetTitle("Multiplicity");
 476   fOutputList->Add(fMultDist3);
 477
 478   TH1F *fMultDist4 = new TH1F("fMultDist4","Multiplicity Distribution",kNbinsM,-.5,kNbinsM-.5);
 479   fMultDist4->GetXaxis()->SetTitle("Multiplicity");
 480   fOutputList->Add(fMultDist4);
 481
 482   TH1F *fMultDist5 = new TH1F("fMultDist5","Multiplicity Distribution",kNbinsM,-.5,kNbinsM-.5);
 483   fMultDist5->GetXaxis()->SetTitle("Multiplicity");
 484   fOutputList->Add(fMultDist5);
 485
 486
 487   TH3F *fPtEtaDist = new TH3F("fPtEtaDist","PtEtaDist",2,-1.1,1.1, 300,0,3., 28,-1.4,1.4);
 488   fOutputList->Add(fPtEtaDist);
 489
 490   TH3F *fPhiPtDist = new TH3F("fPhiPtDist","PhiPtDist",2,-1.1,1.1, 120,0,2*PI, 300,0,3.);
 491   fOutputList->Add(fPhiPtDist);
 492
 493
 494   for(Int_t cv=0; cv<kNCutVariations; cv++){
 495
 496     if(cv==0){
 497       TString *nameXi=new TString("fXi_");
 498       TString *nameXibar=new TString("fXibar_");
 499       *nameXi += cv;
 500       *nameXibar += cv;
 501       CutVar[cv].fXi = new TH3F(nameXi->Data(),"Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.2,1.4);
 502       fOutputList->Add(CutVar[cv].fXi);
 503       CutVar[cv].fXibar = new TH3F(nameXibar->Data(),"Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.2,1.4);
 504       fOutputList->Add(CutVar[cv].fXibar);
 505       //
 506       TString *nameMCrecXi = new TString("fMCrecXi_");
 507       TString *nameMCrecXibar = new TString("fMCrecXi_");
 508       *nameMCrecXi += cv;
 509       *nameMCrecXibar += cv;
 510       CutVar[cv].fMCrecXi = new TH3F(nameMCrecXi->Data(),"Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.2,1.4);
 511       CutVar[cv].fMCrecXibar = new TH3F(nameMCrecXibar->Data(),"Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.2,1.4);
 512       fOutputList->Add(CutVar[cv].fMCrecXi);
 513       fOutputList->Add(CutVar[cv].fMCrecXibar);
 514     }
 515     //
 516     TString *nameXiMinusPiPlus = new TString("fXiMinusPiPlus_");
 517     TString *nameXiMinusPiMinus = new TString("fXiMinusPiMinus_");
 518     TString *nameXiPlusPiPlus = new TString("fXiPlusPiPlus_");
 519     TString *nameXiPlusPiMinus = new TString("fXiPlusPiMinus_");
 520     TString *nameXiMinusPiPlusbkg = new TString("fXiMinusPiPlusbkg_");
 521     TString *nameXiMinusPiMinusbkg = new TString("fXiMinusPiMinusbkg_");
 522     TString *nameXiPlusPiPlusbkg = new TString("fXiPlusPiPlusbkg_");
 523     TString *nameXiPlusPiMinusbkg = new TString("fXiPlusPiMinusbkg_");
 524     *nameXiMinusPiPlus += cv;
 525     *nameXiMinusPiMinus += cv;
 526     *nameXiPlusPiPlus += cv;
 527     *nameXiPlusPiMinus += cv;
 528     *nameXiMinusPiPlusbkg += cv;
 529     *nameXiMinusPiMinusbkg += cv;
 530     *nameXiPlusPiPlusbkg += cv;
 531     *nameXiPlusPiMinusbkg += cv;
 532     CutVar[cv].fXiMinusPiPlus  = new TH3F(nameXiMinusPiPlus->Data(),"Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.4,1.6);
 533     CutVar[cv].fXiMinusPiMinus = new TH3F(nameXiMinusPiMinus->Data(),"Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.4,1.6);
 534     CutVar[cv].fXiPlusPiPlus   = new TH3F(nameXiPlusPiPlus->Data(),"Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.4,1.6);
 535     CutVar[cv].fXiPlusPiMinus  = new TH3F(nameXiPlusPiMinus->Data(),"Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.4,1.6);
 536     CutVar[cv].fXiMinusPiPlusbkg  = new TH3F(nameXiMinusPiPlusbkg->Data(),"Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.4,1.6);
 537     CutVar[cv].fXiMinusPiMinusbkg = new TH3F(nameXiMinusPiMinusbkg->Data(),"Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.4,1.6);
 538     CutVar[cv].fXiPlusPiPlusbkg   = new TH3F(nameXiPlusPiPlusbkg->Data(),"Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.4,1.6);
 539     CutVar[cv].fXiPlusPiMinusbkg  = new TH3F(nameXiPlusPiMinusbkg->Data(),"Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.4,1.6);
 540
 541     fOutputList->Add(CutVar[cv].fXiMinusPiPlus);
 542     fOutputList->Add(CutVar[cv].fXiMinusPiMinus);
 543     fOutputList->Add(CutVar[cv].fXiPlusPiPlus);
 544     fOutputList->Add(CutVar[cv].fXiPlusPiMinus);
 545     fOutputList->Add(CutVar[cv].fXiMinusPiPlusbkg);
 546     fOutputList->Add(CutVar[cv].fXiMinusPiMinusbkg);
 547     fOutputList->Add(CutVar[cv].fXiPlusPiPlusbkg);
 548     fOutputList->Add(CutVar[cv].fXiPlusPiMinusbkg);
 549     //
 550
 551
 552     TString *nameMCrecXiMinusPiPlus = new TString("fMCrecXiMinusPiPlus_");
 553     TString *nameMCrecXiPlusPiMinus = new TString("fMCrecXiPlusPiMinus_");
 554     *nameMCrecXiMinusPiPlus += cv;
 555     *nameMCrecXiPlusPiMinus += cv;
 556     CutVar[cv].fMCrecXiMinusPiPlus  = new TH3F(nameMCrecXiMinusPiPlus->Data(),"Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.4,1.6);
 557     CutVar[cv].fMCrecXiPlusPiMinus  = new TH3F(nameMCrecXiPlusPiMinus->Data(),"Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.4,1.6);
 558     fOutputList->Add(CutVar[cv].fMCrecXiMinusPiPlus);
 559     fOutputList->Add(CutVar[cv].fMCrecXiPlusPiMinus);
 560     //
 561
 562     /*
 563     CutVar[cv].fMCrecXiStarxiy = new TH2F("fMCrecXiStarxiy","y distribution",80,-2,2, 80,-2,2);
 564     CutVar[cv].fMCrecXiStarpiony = new TH2F("fMCrecXiStarpiony","y distribution",80,-2,2, 80,-2,2);
 565     fOutputList->Add(CutVar[cv].fMCrecXiStarxiy);
 566     fOutputList->Add(CutVar[cv].fMCrecXiStarpiony);
 567     CutVar[cv].fMCrecXilambday = new TH2F("fMCrecXilambday","y distribution",80,-2,2, 80,-2,2);
 568     CutVar[cv].fMCrecXipiony = new TH2F("fMCrecXipiony","y distribution",80,-2,2, 80,-2,2);
 569     fOutputList->Add(CutVar[cv].fMCrecXilambday);
 570     fOutputList->Add(CutVar[cv].fMCrecXipiony);
 571     CutVar[cv].fMCrecLamprotony = new TH2F("fMCrecLamprotony","y distribution",80,-2,2, 80,-2,2);
 572     CutVar[cv].fMCrecLampiony = new TH2F("fMCrecLampiony","y distribution",80,-2,2, 80,-2,2);
 573     fOutputList->Add(CutVar[cv].fMCrecLamprotony);
 574     fOutputList->Add(CutVar[cv].fMCrecLampiony);
 575     */
 576   }
 577
 578
 579
 580
 581   //////////////////////
 582   // MC input histos
 583   TH3F *fMCinputXiStar = new TH3F("fMCinputXiStar","Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.4,1.6);
 584   TH3F *fMCinputXiStarbar = new TH3F("fMCinputXiStarbar","Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.4,1.6);
 585   fOutputList->Add(fMCinputXiStar);
 586   fOutputList->Add(fMCinputXiStarbar);
 587
 588   TH3F *fMCinputXi = new TH3F("fMCinputXi","Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.2,1.4);
 589   TH3F *fMCinputXibar = new TH3F("fMCinputXibar","Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.2,1.4);
 590   fOutputList->Add(fMCinputXi);
 591   fOutputList->Add(fMCinputXibar);
 592
 593   //
 594
 595   TH3F *fMCinputTotalXiStar1 = new TH3F("fMCinputTotalXiStar1","Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.4,1.6);
 596   TH3F *fMCinputTotalXiStarbar1 = new TH3F("fMCinputTotalXiStarbar1","Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.4,1.6);
 597   fOutputList->Add(fMCinputTotalXiStar1);
 598   fOutputList->Add(fMCinputTotalXiStarbar1);
 599
 600   TH3F *fMCinputTotalXi1 = new TH3F("fMCinputTotalXi1","Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.2,1.4);
 601   TH3F *fMCinputTotalXibar1 = new TH3F("fMCinputTotalXibar1","Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.2,1.4);
 602   fOutputList->Add(fMCinputTotalXi1);
 603   fOutputList->Add(fMCinputTotalXibar1);
 604
 605   //
 606
 607   TH3F *fMCinputTotalXiStar3 = new TH3F("fMCinputTotalXiStar3","Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.4,1.6);
 608   TH3F *fMCinputTotalXiStarbar3 = new TH3F("fMCinputTotalXiStarbar3","Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.4,1.6);
 609   fOutputList->Add(fMCinputTotalXiStar3);
 610   fOutputList->Add(fMCinputTotalXiStarbar3);
 611
 612   TH3F *fMCinputTotalXi3 = new TH3F("fMCinputTotalXi3","Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.2,1.4);
 613   TH3F *fMCinputTotalXibar3 = new TH3F("fMCinputTotalXibar3","Invariant Mass Distribution", 100,0,10, 40,-2,2, 200,1.2,1.4);
 614   fOutputList->Add(fMCinputTotalXi3);
 615   fOutputList->Add(fMCinputTotalXibar3);
 616
 617   //
 618
 619   TH2F *fMCinputXiStarxiy = new TH2F("fMCinputXiStarxiy","y distribution",80,-2,2, 80,-2,2);
 620   TH2F *fMCinputXiStarpiony = new TH2F("fMCinputXiStarpiony","y distribution",80,-2,2, 80,-2,2);
 621   fOutputList->Add(fMCinputXiStarxiy);
 622   fOutputList->Add(fMCinputXiStarpiony);
 623   TH2F *fMCinputXilambday = new TH2F("fMCinputXilambday","y distribution",80,-2,2, 80,-2,2);
 624   TH2F *fMCinputXipiony = new TH2F("fMCinputXipiony","y distribution",80,-2,2, 80,-2,2);
 625   fOutputList->Add(fMCinputXilambday);
 626   fOutputList->Add(fMCinputXipiony);
 627   TH2F *fMCinputLamprotony = new TH2F("fMCinputLamprotony","y distribution",80,-2,2, 80,-2,2);
 628   TH2F *fMCinputLampiony = new TH2F("fMCinputLampiony","y distribution",80,-2,2, 80,-2,2);
 629   fOutputList->Add(fMCinputLamprotony);
 630   fOutputList->Add(fMCinputLampiony);
 631
 632
 633   ///////////////////////////////////
 634
 635   PostData(1, fOutputList);
 636 }
 637
 638 //________________________________________________________________________
 639 void AliXiStar::Exec(Option_t *)
 640 {
 641   // Main loop
 642   // Called for each event
 643   cout<<"===========  Event # "<<fEventCounter+1<<"  ==========="<<endl;
 644   fEventCounter++;
 645
 646   if(fAODcase) {cout<<"AODs not fully supported! Exiting event."<<endl; return;}
 647
 648   if(fAODcase) fAOD = dynamic_cast<AliAODEvent*> (InputEvent());
 649   else fESD = dynamic_cast<AliESDEvent*> (InputEvent());
 650
 651   if(fAODcase) {if (!fAOD) {Printf("ERROR: fAOD not available"); return;}}
 652   else {if (!fESD) {Printf("ERROR: fESD not available"); return;}}
 653
 654
 655   // ESD Trigger Cut
 656   if(!fAODcase){
 657     if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected())) {cout<<"Event Rejected"<<endl; return;}
 658   }
 659
 660   ///////////////////////////////////////////////////////////
 661   const AliAODVertex *PrimaryVertexAOD;
 662   const AliESDVertex *PrimaryVertexESD;
 663
 664   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 665   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 666   fPIDResponse = inputHandler->GetPIDResponse();
 667
 668
 669   TClonesArray *mcArray       = 0x0;
 670   //AliAODMCParticle *mcXi;
 671   //AliAODMCParticle *mcXiStarD2;
 672   //AliAODMCParticle *mcXiStar;
 673   AliMCEvent  *mcEvent        = 0x0;
 674   AliStack    *mcstack        = 0x0;
 675   TParticle   *MCLamD1esd     = 0x0;
 676   TParticle   *MCLamD2esd     = 0x0;
 677   TParticle   *MCLamesd       = 0x0;
 678   TParticle   *MCXiD2esd      = 0x0;
 679   TParticle   *MCXiesd        = 0x0;
 680   TParticle   *MCXiStarD2esd  = 0x0;
 681   TParticle   *MCXiStaresd    = 0x0;
 682
 683   Double_t px1,py1,pz1, px2,py2,pz2;
 684   Double_t p1sq,p2sq,e1,e2,angle;
 685   Double_t dca3d;
 686   Float_t dca2[2];
 687   Double_t xiVtx[3];//, xiStarVtx[3];
 688   Double_t xiP[3], xiStarP[3];
 689   Double_t xiStarMom;
 690   Double_t xiMass, xiStarMass;
 691   Double_t xiPt, xiStarPt;
 692   Double_t xiY, xiStarY;
 693   Double_t xiCharge;
 694   Double_t decayLengthXY;
 695   Double_t pDaughter1[3];
 696   Double_t pDaughter2[3];
 697   Double_t xDaughter1[3];
 698   Double_t xDaughter2[3];
 699   //
 700   Double_t bField=0;
 701   UInt_t status=0;
 702   Int_t positiveTracks=0, negativeTracks=0;
 703   Int_t myTracks=0;
 704   //
 705   Double_t primaryVtx[3]={0};
 706   Int_t mBin=0;
 707   Int_t zBin=0;
 708   Double_t zStep=2*10/Double_t(fZvertexBins), zStart=-10.;
 709   //
 710   Bool_t mcXiFilled=kFALSE;// So that mctracks are never used uninitialized
 711
 712   if(fMCcase){
 713     if(fAODcase){
 714       mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName());
 715       if(!mcArray){
 716         cout<<"No MC particle branch found"<<endl;
 717         return;
 718       }
 719     }else {
 720       mcEvent = MCEvent();
 721       if (!mcEvent) {Printf("ERROR: Could not retrieve MC event"); return;}
 722
 723       mcstack = mcEvent->Stack();
 724       if (!mcstack) {Printf("ERROR: Could not retrieve the stack"); return;}
 725     }
 726   }
 727
 728
 729   /////////////////////////////////////////////////
 730
 731
 732
 733   if(fAODcase){// AOD case
 734     ////////////////////////////////
 735     // Vertexing
 736     ((TH1F*)fOutputList->FindObject("fMultDist1"))->Fill(fAOD->GetNumberOfTracks());
 737     PrimaryVertexAOD = fAOD->GetPrimaryVertex();
 738     if(!PrimaryVertexAOD) return;
 739
 740     if(fMCcase){
 741       for (Int_t it = 0; it < mcArray->GetEntriesFast(); it++) {
 742         AliAODMCParticle *mcInputTrack = (AliAODMCParticle*)mcArray->At(it);
 743         if (!mcInputTrack) {
 744           Error("UserExec", "Could not receive track %d", it);
 745           continue;
 746         }
 747
 748         if(!mcInputTrack->IsPhysicalPrimary()) continue;
 749
 750         // Xi
 751         if(mcInputTrack->GetPdgCode() == +kXiCode) ((TH3F*)fOutputList->FindObject("fMCinputTotalXi1"))->Fill(mcInputTrack->Pt(), mcInputTrack->Y(), mcInputTrack->GetCalcMass());
 752         if(mcInputTrack->GetPdgCode() == -kXiCode) ((TH3F*)fOutputList->FindObject("fMCinputTotalXibar1"))->Fill(mcInputTrack->Pt(), mcInputTrack->Y(), mcInputTrack->GetCalcMass());
 753
 754         // XiStar
 755         if(mcInputTrack->GetPdgCode() == +kXiStarCode) ((TH3F*)fOutputList->FindObject("fMCinputTotalXiStar1"))->Fill(mcInputTrack->Pt(), mcInputTrack->Y(), mcInputTrack->GetCalcMass());
 756         if(mcInputTrack->GetPdgCode() == -kXiStarCode) ((TH3F*)fOutputList->FindObject("fMCinputTotalXiStarbar1"))->Fill(mcInputTrack->Pt(), mcInputTrack->Y(), mcInputTrack->GetCalcMass());
 757       }
 758     }
 759
 760
 761
 762     primaryVtx[0]=PrimaryVertexAOD->GetX(); primaryVtx[1]=PrimaryVertexAOD->GetY(); primaryVtx[2]=PrimaryVertexAOD->GetZ();
 763     ((TH3F*)fOutputList->FindObject("fVertexDist1"))->Fill(primaryVtx[0], primaryVtx[1], primaryVtx[2]);
 764     if(fabs(primaryVtx[2]) > 10.) return; // Z-Vertex Cut
 765     ((TH1F*)fOutputList->FindObject("fMultDist2"))->Fill(fAOD->GetNumberOfTracks());
 766
 767
 768     if(fAOD->IsPileupFromSPD()) return; // Reject Pile-up events
 769
 770     ((TH1F*)fOutputList->FindObject("fMultDist3"))->Fill(fAOD->GetNumberOfTracks());
 771     ((TH3F*)fOutputList->FindObject("fVertexDist3"))->Fill(primaryVtx[0], primaryVtx[1], primaryVtx[2]);
 772
 773     if(PrimaryVertexAOD->GetNContributors() >= 1) ((TH1F*)fOutputList->FindObject("fMultDist4"))->Fill(fAOD->GetNumberOfTracks());
 774
 775
 776
 777     Printf("There are %d tracks in this event", fAOD->GetNumberOfTracks());
 778     // fNtracks Cut
 779     if(fAOD->GetNumberOfTracks() > kNbinsM) {cout<<"More tracks than limit"<<endl; return;}
 780     bField = fAOD->GetMagneticField();
 781
 782     // Track loop
 783     for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) {
 784       AliAODTrack* aodtrack = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(i));
 785       if(!aodtrack) AliFatal("Not a standard AOD");
 786       if (!aodtrack) continue;
 787
 788       status=aodtrack->GetStatus();
 789
 790
 791       if( (status&AliESDtrack::kTPCrefit)==0) continue;// Require tpcrefit
 792       if( aodtrack->GetTPCNcls() < 70) continue;// TPC Ncluster cut
 793       if(aodtrack->Pt() < 0.15) continue;
 794       AliAODVertex *VtxType=aodtrack->GetProdVertex();
 795       if((VtxType->GetType())==AliAODVertex::kKink) continue;// Reject Kinks
 796
 797       Bool_t goodMomentum = aodtrack->GetPxPyPz( fTempStruct[myTracks].fP);
 798       if(!goodMomentum) continue;
 799       aodtrack->GetXYZ( fTempStruct[myTracks].fX);
 800
 801
 802       aodtrack->GetCovarianceXYZPxPyPz( fTempStruct[myTracks].fCov);
 803
 804       dca2[0] = sqrt( pow(fTempStruct[myTracks].fX[0] - primaryVtx[0],2) + pow(fTempStruct[myTracks].fX[1] - primaryVtx[1],2));
 805       dca2[1] = sqrt( pow(fTempStruct[myTracks].fX[2] - primaryVtx[2],2));
 806       dca3d = sqrt( pow(dca2[0],2) + pow(dca2[1],2));
 807
 808       ((TH1F*)fOutputList->FindObject("fDCADist"))->Fill(fAOD->GetNumberOfTracks(), dca3d);
 809       ((TH3F*)fOutputList->FindObject("fPhiPtDist"))->Fill(aodtrack->Charge(), aodtrack->Phi(), aodtrack->Pt());
 810       ((TH3F*)fOutputList->FindObject("fPtEtaDist"))->Fill(aodtrack->Charge(), aodtrack->Pt(), aodtrack->Eta());
 811
 812
 813
 814
 815       fTempStruct[myTracks].fStatus = status;
 816       fTempStruct[myTracks].fFilterMap = aodtrack->GetFilterMap();
 817       fTempStruct[myTracks].fID = aodtrack->GetID();
 818       fTempStruct[myTracks].fLabel = aodtrack->GetLabel();
 819       fTempStruct[myTracks].fPhi = atan2(fTempStruct[myTracks].fP[1], fTempStruct[myTracks].fP[0]);
 820       if(fTempStruct[myTracks].fPhi < 0) fTempStruct[myTracks].fPhi += 2*PI;
 821       fTempStruct[myTracks].fPt = sqrt(pow(fTempStruct[myTracks].fP[0],2) + pow(fTempStruct[myTracks].fP[1],2));
 822       fTempStruct[myTracks].fMom = sqrt( pow(fTempStruct[myTracks].fPt,2) + pow(fTempStruct[myTracks].fP[2],2) );
 823       fTempStruct[myTracks].fEta = aodtrack->Eta();
 824       fTempStruct[myTracks].fCharge = aodtrack->Charge();
 825       fTempStruct[myTracks].fDCAXY = dca2[0];
 826       fTempStruct[myTracks].fDCAZ = dca2[1];
 827       fTempStruct[myTracks].fDCA = dca3d;
 828       fTempStruct[myTracks].fNSigmaPi = fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kPion));
 829       fTempStruct[myTracks].fNSigmaK = fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kKaon));
 830       fTempStruct[myTracks].fNSigmaPr = fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kProton));
 831
 832
 833       if(aodtrack->Charge() > 0) positiveTracks++;
 834       else negativeTracks++;
 835
 836
 837       myTracks++;
 838     }
 839   }else {// ESDs
 840
 841     ((TH1F*)fOutputList->FindObject("fMultDist1"))->Fill(fESD->GetNumberOfTracks());
 842     PrimaryVertexESD = fESD->GetPrimaryVertex();
 843     if(!PrimaryVertexESD) return;
 844
 845     primaryVtx[0]=PrimaryVertexESD->GetX(); primaryVtx[1]=PrimaryVertexESD->GetY(); primaryVtx[2]=PrimaryVertexESD->GetZ();
 846     ((TH3F*)fOutputList->FindObject("fVertexDist1"))->Fill(primaryVtx[0], primaryVtx[1], primaryVtx[2]);
 847
 848     if(fMCcase){
 849       /////////////////////////////////////////////////
 850       // Lam mc input
 851       /////////////////////////////////////////////////
 852       for (Int_t it = 0; it < mcstack->GetNprimary(); it++) {
 853         TParticle *mcInputTrack = (TParticle*)mcstack->Particle(it);
 854         if (!mcInputTrack) {
 855           Error("UserExec", "Could not receive track %d", it);
 856           continue;
 857         }
 858
 859
 860         // Xi
 861         if(mcInputTrack->GetPdgCode() == +kXiCode) ((TH3F*)fOutputList->FindObject("fMCinputTotalXi1"))->Fill(mcInputTrack->Pt(), mcInputTrack->Y(), mcInputTrack->GetCalcMass());
 862         if(mcInputTrack->GetPdgCode() == -kXiCode) ((TH3F*)fOutputList->FindObject("fMCinputTotalXibar1"))->Fill(mcInputTrack->Pt(), mcInputTrack->Y(), mcInputTrack->GetCalcMass());
 863
 864         // XiStar
 865         if(mcInputTrack->GetPdgCode() == +kXiStarCode) ((TH3F*)fOutputList->FindObject("fMCinputTotalXiStar1"))->Fill(mcInputTrack->Pt(), mcInputTrack->Y(), mcInputTrack->GetCalcMass());
 866         if(mcInputTrack->GetPdgCode() == -kXiStarCode) ((TH3F*)fOutputList->FindObject("fMCinputTotalXiStarbar1"))->Fill(mcInputTrack->Pt(), mcInputTrack->Y(), mcInputTrack->GetCalcMass());
 867
 868       }
 869
 870
 871     }
 872
 873
 874
 875     if(fabs(primaryVtx[2]) > 10.) return; // Z-Vertex Cut
 876     ((TH1F*)fOutputList->FindObject("fMultDist2"))->Fill(fESD->GetNumberOfTracks());
 877
 878
 879     if(fESD->IsPileupFromSPD()) return; // Reject Pile-up events
 880
 881     ((TH1F*)fOutputList->FindObject("fMultDist3"))->Fill(fESD->GetNumberOfTracks());
 882     ((TH3F*)fOutputList->FindObject("fVertexDist3"))->Fill(primaryVtx[0], primaryVtx[1], primaryVtx[2]);
 883
 884     if(PrimaryVertexESD->GetNContributors() >= 1) ((TH1F*)fOutputList->FindObject("fMultDist4"))->Fill(fESD->GetNumberOfTracks());
 885
 886     Printf("There are %d tracks in this event", fESD->GetNumberOfTracks());
 887     // fNtracks Cut
 888     if(fESD->GetNumberOfTracks() > kNbinsM) {cout<<"More tracks than limit"<<endl; return;}
 889     bField = fESD->GetMagneticField();
 890
 891
 892     // Track loop
 893     for (Int_t i = 0; i < fESD->GetNumberOfTracks(); i++) {
 894       AliESDtrack* esdtrack = fESD->GetTrack(i);
 895       if (!esdtrack) continue;
 896
 897       status=esdtrack->GetStatus();
 898
 899       if(!fTrackCut->AcceptTrack(esdtrack)) continue;
 900
 901       Bool_t goodMomentum = esdtrack->GetPxPyPz( fTempStruct[myTracks].fP);
 902       if(!goodMomentum) continue;
 903       esdtrack->GetXYZ( fTempStruct[myTracks].fX);
 904
 905
 906       esdtrack->GetCovarianceXYZPxPyPz( fTempStruct[myTracks].fCov);
 907       //esdtrack->GetImpactParameters(dca2, cov);
 908       dca2[0] = sqrt( pow(fTempStruct[myTracks].fX[0] - primaryVtx[0],2) + pow(fTempStruct[myTracks].fX[1] - primaryVtx[1],2));
 909       dca2[1] = sqrt( pow(fTempStruct[myTracks].fX[2] - primaryVtx[2],2));
 910       dca3d = sqrt( pow(dca2[0],2) + pow(dca2[1],2));
 911
 912       ((TH1F*)fOutputList->FindObject("fDCADist"))->Fill(fESD->GetNumberOfTracks(), dca3d);
 913       ((TH3F*)fOutputList->FindObject("fPhiPtDist"))->Fill(esdtrack->Charge(), esdtrack->Phi(), esdtrack->Pt());
 914       ((TH3F*)fOutputList->FindObject("fPtEtaDist"))->Fill(esdtrack->Charge(), esdtrack->Pt(), esdtrack->Eta());
 915
 916
 917
 918       fTempStruct[myTracks].fStatus = status;
 919       fTempStruct[myTracks].fID = esdtrack->GetID();
 920       fTempStruct[myTracks].fLabel = esdtrack->GetLabel();
 921       fTempStruct[myTracks].fPhi = atan2(fTempStruct[myTracks].fP[1], fTempStruct[myTracks].fP[0]);
 922       if(fTempStruct[myTracks].fPhi < 0) fTempStruct[myTracks].fPhi += 2*PI;
 923       fTempStruct[myTracks].fPt = sqrt(pow(fTempStruct[myTracks].fP[0],2) + pow(fTempStruct[myTracks].fP[1],2));
 924       fTempStruct[myTracks].fMom = sqrt( pow(fTempStruct[myTracks].fPt,2) + pow(fTempStruct[myTracks].fP[2],2) );
 925       fTempStruct[myTracks].fEta = esdtrack->Eta();
 926       fTempStruct[myTracks].fCharge = esdtrack->Charge();
 927       fTempStruct[myTracks].fDCAXY = dca2[0];
 928       fTempStruct[myTracks].fDCAZ = dca2[1];
 929       fTempStruct[myTracks].fDCA = dca3d;
 930       fTempStruct[myTracks].fNSigmaPi = fabs(fPIDResponse->NumberOfSigmasTPC(esdtrack,AliPID::kPion));
 931       fTempStruct[myTracks].fNSigmaK = fabs(fPIDResponse->NumberOfSigmasTPC(esdtrack,AliPID::kKaon));
 932       fTempStruct[myTracks].fNSigmaPr = fabs(fPIDResponse->NumberOfSigmasTPC(esdtrack,AliPID::kProton));
 933       fTempStruct[myTracks].fNclusTPC = esdtrack->GetTPCNcls();
 934
 935
 936       if(esdtrack->Charge() > 0) positiveTracks++;
 937       else negativeTracks++;
 938
 939       myTracks++;
 940     }
 941
 942   }// end of ESD case
 943
 944
 945
 946   if(myTracks >= 1) {
 947     ((TH1F*)fOutputList->FindObject("fMultDist5"))->Fill(myTracks);
 948     ((TH3F*)fOutputList->FindObject("fMultDist3d"))->Fill(positiveTracks+negativeTracks, positiveTracks, negativeTracks);
 949   }
 950
 951
 952   cout<<"There are "<<myTracks<<"  myTracks"<<endl;
 953
 954   // set Z Vertex bin
 955   for(Int_t i=0; i<fZvertexBins; i++){
 956     if( (primaryVtx[2] > zStart+i*zStep) && (primaryVtx[2] < zStart+(i+1)*zStep) ){
 957       zBin=i;
 958       break;
 959     }
 960   }
 961
 962   // set Multiplicity bin
 963   for(Int_t i=0; i<fMultBins; i++){
 964     if( ( myTracks > fMultLimits[i]) && ( myTracks <= fMultLimits[i+1]) ) { mBin=i; break;}
 965   }
 966
 967
 968
 969   ////////////////////////////////////
 970   // Add event to buffer if > 0 tracks
 971   if(myTracks > 0){
 972     fEC[zBin][mBin]->FIFOShift();
 973     (fEvt) = fEC[zBin][mBin]->fEvtStr;
 974     (fEvt)->fNTracks = myTracks;
 975     for(Int_t i=0; i<myTracks; i++) (fEvt)->fTracks[i] = fTempStruct[i];
 976   }
 977
 978
 979
 980
 981   if(fMCcase && fAODcase){// get Input MC information for AOD case
 982
 983     /////////////////////////////////////////////////
 984     // Xi mc input
 985     /////////////////////////////////////////////////
 986     for (Int_t it = 0; it < mcArray->GetEntriesFast(); it++) {
 987       AliAODMCParticle *mcInputTrackXi = (AliAODMCParticle*)mcArray->At(it);
 988
 989       if (!mcInputTrackXi) {
 990         Error("UserExec", "Could not receive track %d", it);
 991         continue;
 992       }
 993
 994       if(abs(mcInputTrackXi->GetPdgCode())!=kXiCode) continue;
 995       if(!mcInputTrackXi->IsPhysicalPrimary()) continue;
 996
 997
 998       if(mcInputTrackXi->GetPdgCode() == +kXiCode) ((TH3F*)fOutputList->FindObject("fMCinputTotalXi3"))->Fill(mcInputTrackXi->Pt(), mcInputTrackXi->Y(), mcInputTrackXi->GetCalcMass());
 999       else ((TH3F*)fOutputList->FindObject("fMCinputTotalXibar3"))->Fill(mcInputTrackXi->Pt(), mcInputTrackXi->Y(), mcInputTrackXi->GetCalcMass());
1000
1001
1002
1003       AliAODMCParticle *mcInputTrackXiD1 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXi->GetDaughter(0)));
1004       AliAODMCParticle *mcInputTrackXiD2 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXi->GetDaughter(1)));
1005
1006       if(abs(mcInputTrackXiD1->GetPdgCode())!=kLambdaCode && abs(mcInputTrackXiD2->GetPdgCode())!=kLambdaCode) continue;
1007       if(abs(mcInputTrackXiD1->GetPdgCode())!=kPionCode && abs(mcInputTrackXiD2->GetPdgCode())!=kPionCode) continue;
1008
1009
1010       // At this point we have the right Xi decay channel
1011
1012       AliAODMCParticle *mcInputTrackLamD1;
1013       AliAODMCParticle *mcInputTrackLamD2;
1014
1015       if(abs(mcInputTrackXiD1->GetPdgCode())==kLambdaCode) {
1016         mcInputTrackLamD1 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiD1->GetDaughter(0)));
1017         mcInputTrackLamD2 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiD1->GetDaughter(1)));
1018       }
1019       else{
1020         mcInputTrackLamD1 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiD2->GetDaughter(0)));
1021         mcInputTrackLamD2 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiD2->GetDaughter(1)));
1022       }
1023
1024
1025       if(abs(mcInputTrackLamD1->GetPdgCode())!=kProtonCode && abs(mcInputTrackLamD2->GetPdgCode())!=kProtonCode) continue;
1026       if(abs(mcInputTrackLamD1->GetPdgCode())!=kPionCode && abs(mcInputTrackLamD2->GetPdgCode())!=kPionCode) continue;
1027
1028       // At this point we have the right Lambda decay channel
1029
1030       if(mcInputTrackXi->GetPdgCode() == +kXiCode) ((TH3F*)fOutputList->FindObject("fMCinputXi"))->Fill(mcInputTrackXi->Pt(), mcInputTrackXi->Y(), mcInputTrackXi->GetCalcMass());
1031       else ((TH3F*)fOutputList->FindObject("fMCinputXibar"))->Fill(mcInputTrackXi->Pt(), mcInputTrackXi->Y(), mcInputTrackXi->GetCalcMass());
1032
1033     }
1034
1035
1036
1037     /////////////////////////////////////////////////
1038     // XiStar mc input
1039     /////////////////////////////////////////////////
1040     for (Int_t it = 0; it < mcArray->GetEntriesFast(); it++) {
1041       AliAODMCParticle *mcInputTrackXiStar = (AliAODMCParticle*)mcArray->At(it);
1042       if (!mcInputTrackXiStar) {
1043         Error("UserExec", "Could not receive track %d", it);
1044         continue;
1045       }
1046
1047       if(abs(mcInputTrackXiStar->GetPdgCode())!=kXiStarCode) continue;
1048       if(!mcInputTrackXiStar->IsPhysicalPrimary()) continue;
1049
1050       if(mcInputTrackXiStar->GetPdgCode() == +kXiStarCode) ((TH3F*)fOutputList->FindObject("fMCinputTotalXiStar3"))->Fill(mcInputTrackXiStar->Pt(), mcInputTrackXiStar->Y(), mcInputTrackXiStar->GetCalcMass());
1051       else ((TH3F*)fOutputList->FindObject("fMCinputTotalXiStarbar3"))->Fill(mcInputTrackXiStar->Pt(), mcInputTrackXiStar->Y(), mcInputTrackXiStar->GetCalcMass());
1052
1053
1054       AliAODMCParticle *mcInputTrackXiStarD1 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiStar->GetDaughter(0)));
1055       AliAODMCParticle *mcInputTrackXiStarD2 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiStar->GetDaughter(1)));
1056
1057       if(abs(mcInputTrackXiStarD1->GetPdgCode())!=kXiCode && abs(mcInputTrackXiStarD2->GetPdgCode())!=kXiCode) continue;
1058       if(abs(mcInputTrackXiStarD1->GetPdgCode())!=kPionCode && abs(mcInputTrackXiStarD2->GetPdgCode())!=kPionCode) continue;
1059
1060
1061       // At this point we have the right Xi(1530) decay channel
1062
1063       AliAODMCParticle *mcInputTrackXiD1;
1064       AliAODMCParticle *mcInputTrackXiD2;
1065       if(abs(mcInputTrackXiStarD1->GetPdgCode())==kXiCode) {
1066         mcInputTrackXiD1 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiStarD1->GetDaughter(0)));
1067         mcInputTrackXiD2 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiStarD1->GetDaughter(1)));
1068       }
1069       else{
1070         mcInputTrackXiD1 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiStarD2->GetDaughter(0)));
1071         mcInputTrackXiD2 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiStarD2->GetDaughter(1)));
1072       }
1073
1074       if(abs(mcInputTrackXiD1->GetPdgCode())!=kLambdaCode && abs(mcInputTrackXiD2->GetPdgCode())!=kLambdaCode) continue;
1075       if(abs(mcInputTrackXiD1->GetPdgCode())!=kPionCode && abs(mcInputTrackXiD2->GetPdgCode())!=kPionCode) continue;
1076
1077
1078       // At this point we have the right Xi decay channel
1079
1080       AliAODMCParticle *mcInputTrackLamD1;
1081       AliAODMCParticle *mcInputTrackLamD2;
1082
1083       if(abs(mcInputTrackXiD1->GetPdgCode())==kLambdaCode) {
1084         mcInputTrackLamD1 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiD1->GetDaughter(0)));
1085         mcInputTrackLamD2 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiD1->GetDaughter(1)));
1086       }
1087       else{
1088         mcInputTrackLamD1 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiD2->GetDaughter(0)));
1089         mcInputTrackLamD2 = (AliAODMCParticle*)mcArray->At(abs(mcInputTrackXiD2->GetDaughter(1)));
1090       }
1091
1092       if(abs(mcInputTrackLamD1->GetPdgCode())!=kProtonCode && abs(mcInputTrackLamD2->GetPdgCode())!=kProtonCode) continue;
1093       if(abs(mcInputTrackLamD1->GetPdgCode())!=kPionCode && abs(mcInputTrackLamD2->GetPdgCode())!=kPionCode) continue;
1094
1095       // At this point we the right Lambda decay channel
1096
1097       if(mcInputTrackXiStar->GetPdgCode() == +kXiStarCode) ((TH3F*)fOutputList->FindObject("fMCinputXiStar"))->Fill(mcInputTrackXiStar->Pt(), mcInputTrackXiStar->Y(), mcInputTrackXiStar->GetCalcMass());
1098       else ((TH3F*)fOutputList->FindObject("fMCinputXiStarbar"))->Fill(mcInputTrackXiStar->Pt(), mcInputTrackXiStar->Y(), mcInputTrackXiStar->GetCalcMass());
1099
1100       if(abs(mcInputTrackXiStarD1->GetPdgCode())==kXiCode) {
1101         ((TH2F*)fOutputList->FindObject("fMCinputXiStarxiy"))->Fill(mcInputTrackXiStar->Y(), mcInputTrackXiStarD1->Y());
1102         ((TH2F*)fOutputList->FindObject("fMCinputXiStarpiony"))->Fill(mcInputTrackXiStar->Y(), mcInputTrackXiStarD2->Y());
1103       }else{
1104         ((TH2F*)fOutputList->FindObject("fMCinputXiStarxiy"))->Fill(mcInputTrackXiStar->Y(), mcInputTrackXiStarD2->Y());
1105         ((TH2F*)fOutputList->FindObject("fMCinputXiStarpiony"))->Fill(mcInputTrackXiStar->Y(), mcInputTrackXiStarD1->Y());
1106       }
1107       if(abs(mcInputTrackXiD1->GetPdgCode())==kLambdaCode){
1108         ((TH2F*)fOutputList->FindObject("fMCinputXilambday"))->Fill(mcInputTrackXiStar->Y(), mcInputTrackXiD1->Y());
1109         ((TH2F*)fOutputList->FindObject("fMCinputXipiony"))->Fill(mcInputTrackXiStar->Y(), mcInputTrackXiD2->Y());
1110       }else{
1111         ((TH2F*)fOutputList->FindObject("fMCinputXilambday"))->Fill(mcInputTrackXiStar->Y(), mcInputTrackXiD2->Y());
1112         ((TH2F*)fOutputList->FindObject("fMCinputXipiony"))->Fill(mcInputTrackXiStar->Y(), mcInputTrackXiD1->Y());
1113       }
1114       if(abs(mcInputTrackLamD1->GetPdgCode())==kProtonCode){
1115         ((TH2F*)fOutputList->FindObject("fMCinputLamprotony"))->Fill(mcInputTrackXiStar->Y(), mcInputTrackLamD1->Y());
1116         ((TH2F*)fOutputList->FindObject("fMCinputLampiony"))->Fill(mcInputTrackXiStar->Y(), mcInputTrackLamD2->Y());
1117       }else {
1118         ((TH2F*)fOutputList->FindObject("fMCinputLamprotony"))->Fill(mcInputTrackXiStar->Y(), mcInputTrackLamD2->Y());
1119         ((TH2F*)fOutputList->FindObject("fMCinputLampiony"))->Fill(mcInputTrackXiStar->Y(), mcInputTrackLamD1->Y());
1120       }
1121
1122
1123     }
1124   }
1125
1126
1127
1128   if(fMCcase && !fAODcase){// get Input MC information for ESD case
1129
1130     /////////////////////////////////////////////////
1131     // Xi mc input
1132     /////////////////////////////////////////////////
1133     for (Int_t it = 0; it < mcstack->GetNprimary(); it++) {
1134       TParticle *mcInputTrackXi = (TParticle*)mcstack->Particle(it);
1135       if (!mcInputTrackXi) {
1136         Error("UserExec", "Could not receive track %d", it);
1137         continue;
1138       }
1139
1140       //if(!mcstack->IsPhysicalPrimary(it)) continue;
1141       if(abs(mcInputTrackXi->GetPdgCode())!=kXiCode) continue;
1142
1143
1144       if(mcInputTrackXi->GetPdgCode() == +kXiCode) ((TH3F*)fOutputList->FindObject("fMCinputTotalXi3"))->Fill(mcInputTrackXi->Pt(), mcInputTrackXi->Y(), mcInputTrackXi->GetCalcMass());
1145       else ((TH3F*)fOutputList->FindObject("fMCinputTotalXibar3"))->Fill(mcInputTrackXi->Pt(), mcInputTrackXi->Y(), mcInputTrackXi->GetCalcMass());
1146
1147     }
1148
1149
1150     /////////////////////////////////////////////////
1151     // XiStar mc input
1152     /////////////////////////////////////////////////
1153     for (Int_t it = 0; it < mcstack->GetNprimary(); it++) {
1154       TParticle *mcInputTrackXiStar = (TParticle*)mcstack->Particle(it);
1155       if (!mcInputTrackXiStar) {
1156         Error("UserExec", "Could not receive track %d", it);
1157         continue;
1158       }
1159
1160       //if(!mcstack->IsPhysicalPrimary(it)) continue;
1161       if(abs(mcInputTrackXiStar->GetPdgCode())!=kXiStarCode) continue;
1162
1163       if(mcInputTrackXiStar->GetPdgCode() == +kXiStarCode) ((TH3F*)fOutputList->FindObject("fMCinputTotalXiStar3"))->Fill(mcInputTrackXiStar->Pt(), mcInputTrackXiStar->Y(), mcInputTrackXiStar->GetCalcMass());
1164       else ((TH3F*)fOutputList->FindObject("fMCinputTotalXiStarbar3"))->Fill(mcInputTrackXiStar->Pt(), mcInputTrackXiStar->Y(), mcInputTrackXiStar->GetCalcMass());
1165
1166     }
1167   }
1168
1169
1170
1171   if(fAODcase) {cout<<"AOD XiVertexer not currently supported! Exiting event"<<endl; return;}
1172
1173   ////////////////////////////////////////////////
1174   // Reconstruction
1175
1176   for(Int_t i=0; i<fESD->GetNumberOfCascades(); i++){
1177
1178     AliESDcascade *Xicandidate = fESD->GetCascade(i);
1179
1180     if(TMath::Abs( Xicandidate->GetPindex()) == TMath::Abs( Xicandidate->GetNindex())) continue;
1181     if(TMath::Abs( Xicandidate->GetPindex()) == TMath::Abs( Xicandidate->GetBindex())) continue;
1182     if(TMath::Abs( Xicandidate->GetNindex()) == TMath::Abs( Xicandidate->GetBindex())) continue;
1183
1184     AliESDtrack *pTrackXi       = fESD->GetTrack(TMath::Abs( Xicandidate->GetPindex()));
1185     AliESDtrack *nTrackXi       = fESD->GetTrack(TMath::Abs( Xicandidate->GetNindex()));
1186     AliESDtrack *bTrackXi       = fESD->GetTrack(TMath::Abs( Xicandidate->GetBindex()));
1187
1188     // Standard track QA cuts
1189     if(!fTrackCut->AcceptTrack(pTrackXi)) continue;
1190     if(!fTrackCut->AcceptTrack(nTrackXi)) continue;
1191     if(!fTrackCut->AcceptTrack(bTrackXi)) continue;
1192
1193     //////////////////////
1194     // DecayParameters Key (number represents array index)
1195     // NclustersTPC: 0=proton, 1=pion first, 2=pion second, 3=pion third
1196     // DCAVtx: 4=proton, 5=pion first, 6=pion second, 7=lambda, 8=pion third
1197     // 9 = DCA proton-pion
1198     // 10 = DCA Lambda-pion
1199     // 11 = Rxy Lambda
1200     // 12 = Rxy Xi
1201     // 13 = Cos PA Lambda
1202     // 14 = Cos PA Xi
1203
1204
1205     fDecayParameters[2] = bTrackXi->GetTPCNcls();
1206
1207
1208     if(Xicandidate->Charge() == -1){
1209       fDecayParameters[0] = pTrackXi->GetTPCNcls();
1210       fDecayParameters[1] = nTrackXi->GetTPCNcls();
1211       fDecayParameters[4] = fabs(pTrackXi->GetD(primaryVtx[0],primaryVtx[1],bField));// DCA Vtx proton
1212       fDecayParameters[5] = fabs(nTrackXi->GetD(primaryVtx[0],primaryVtx[1],bField));// DCA Vtx pion first
1213     }else{
1214       fDecayParameters[0] = nTrackXi->GetTPCNcls();
1215       fDecayParameters[1] = pTrackXi->GetTPCNcls();
1216       fDecayParameters[4] = fabs(nTrackXi->GetD(primaryVtx[0],primaryVtx[1],bField));// DCA Vtx proton
1217       fDecayParameters[5] = fabs(pTrackXi->GetD(primaryVtx[0],primaryVtx[1],bField));// DCA Vtx pion first
1218     }
1219
1220
1221     fDecayParameters[6] = fabs(bTrackXi->GetD(primaryVtx[0],primaryVtx[1],bField));// DCA Vtx pion second
1222     fDecayParameters[7] = fabs(Xicandidate->GetD(primaryVtx[0],primaryVtx[1],primaryVtx[2]));// DCA Vtx Lambda
1223     fDecayParameters[9] = fabs(Xicandidate->GetDcaV0Daughters());// DCA proton-pion
1224     fDecayParameters[10] = fabs(Xicandidate->GetDcaXiDaughters());// DCA Lambda-pion
1225
1226
1227
1228     Double_t tempX[3]={0};
1229     Xicandidate->GetXYZ(tempX[0], tempX[1], tempX[2]);
1230     fDecayParameters[11] = sqrt( pow(tempX[0],2) + pow(tempX[1],2));// Rxy Lambda
1231     if(sqrt( pow(tempX[0],2) + pow(tempX[1],2) ) > fMaxDecayLength) continue;
1232
1233
1234     fDecayParameters[13] = Xicandidate->GetV0CosineOfPointingAngle(primaryVtx[0],primaryVtx[1],primaryVtx[2]);// Cos PA Lambda
1235     fDecayParameters[14] = Xicandidate->GetCascadeCosineOfPointingAngle(primaryVtx[0],primaryVtx[1],primaryVtx[2]);// Cos PA Xi
1236
1237     xiP[0] = Xicandidate->Px();
1238     xiP[1] = Xicandidate->Py();
1239     xiP[2] = Xicandidate->Pz();
1240     xiVtx[0] = Xicandidate->Xv();
1241     xiVtx[1] = Xicandidate->Yv();
1242     xiVtx[2] = Xicandidate->Zv();
1243     xiPt = Xicandidate->Pt();
1244     xiY = Xicandidate->RapXi();
1245     xiMass = Xicandidate->M();
1246     xiCharge = Xicandidate->Charge();
1247
1248     decayLengthXY = sqrt( pow(xiVtx[0]-primaryVtx[0],2) + pow(xiVtx[1]-primaryVtx[1],2) );
1249     fDecayParameters[12] = decayLengthXY;// Rxy Xi
1250     if(decayLengthXY > fMaxDecayLength) continue;// 2d version
1251
1252     Bool_t StandardXi=kTRUE;
1253     if(fDecayParameters[0] < fCutValues[0][0]) StandardXi=kFALSE;// Nclus proton
1254     if(fDecayParameters[1] < fCutValues[0][1]) StandardXi=kFALSE;// Nclus pion first
1255     if(fDecayParameters[2] < fCutValues[0][2]) StandardXi=kFALSE;// Nclus pion second
1256     //
1257     if(fDecayParameters[4] < fCutValues[0][4]) StandardXi=kFALSE;// DCAVtx proton
1258     if(fDecayParameters[5] < fCutValues[0][5]) StandardXi=kFALSE;// DCAVtx pion first
1259     if(fDecayParameters[6] < fCutValues[0][6]) StandardXi=kFALSE;// DCAVtx pion second
1260     if(fDecayParameters[7] < fCutValues[0][7]) StandardXi=kFALSE;// DCAVtx Lambda
1261     //
1262     if(fDecayParameters[9] > fCutValues[0][9]) StandardXi=kFALSE;// DCAV proton-pion
1263     if(fDecayParameters[10] > fCutValues[0][10]) StandardXi=kFALSE;// DCAV Lambda-pion
1264     //
1265     if(fDecayParameters[11] < fCutValues[0][11]) StandardXi=kFALSE;// Rxy Lambda
1266     if(fDecayParameters[12] < fCutValues[0][12]) StandardXi=kFALSE;// Rxy Xi
1267     //
1268     if(fDecayParameters[13] < fCutValues[0][13]) StandardXi=kFALSE;// Cos PA Lambda
1269     if(fDecayParameters[14] < fCutValues[0][14]) StandardXi=kFALSE;// Cos PA Xi
1270
1271     if(StandardXi){
1272       if(xiCharge == -1) CutVar[0].fXi->Fill(xiPt, xiY, xiMass);
1273       else CutVar[0].fXibar->Fill(xiPt, xiY, xiMass);
1274     }
1275
1276     // MC associaton
1277     mcXiFilled = kFALSE;
1278     if(fMCcase && !fAODcase){
1279
1280       MCXiD2esd = (TParticle*)mcstack->Particle(abs(bTrackXi->GetLabel()));
1281
1282       if(abs(MCXiD2esd->GetPdgCode())==kPionCode){
1283
1284         MCLamD1esd = (TParticle*)mcstack->Particle(abs(pTrackXi->GetLabel()));
1285         MCLamD2esd = (TParticle*)mcstack->Particle(abs(nTrackXi->GetLabel()));
1286
1287         if(MCLamD1esd->GetMother(0) == MCLamD2esd->GetMother(0)){
1288           if(abs(MCLamD1esd->GetPdgCode())==kProtonCode || abs(MCLamD2esd->GetPdgCode())==kProtonCode) {
1289             if(abs(MCLamD1esd->GetPdgCode())==kPionCode || abs(MCLamD2esd->GetPdgCode())==kPionCode) {
1290
1291               MCLamesd = (TParticle*)mcstack->Particle(abs(MCLamD1esd->GetMother(0)));
1292               if(abs(MCLamesd->GetPdgCode())==kLambdaCode) {
1293
1294                 if(MCLamesd->GetMother(0) == MCXiD2esd->GetMother(0)){
1295                   MCXiesd = (TParticle*)mcstack->Particle(abs(MCLamesd->GetMother(0)));
1296                   if(abs(MCXiesd->GetPdgCode())==kXiCode) {
1297                     mcXiFilled = kTRUE;
1298
1299                     if(StandardXi){
1300                       if(Xicandidate->Charge() == -1) {
1301                         CutVar[0].fMCrecXi->Fill(xiPt, xiY, xiMass);
1302                       }else {
1303                         CutVar[0].fMCrecXibar->Fill(xiPt, xiY, xiMass);
1304                       }
1305                     }
1306
1307                   }
1308                 }
1309               }
1310             }
1311           }
1312         }
1313       }
1314     }// MC association
1315
1316
1317     if(fabs(xiMass-fTrueMassXi) > fMassWindow) continue;
1318
1319
1320
1321     fXiTrack->Set(xiVtx, xiP, fCovMatrix, Short_t(xiCharge));
1322
1323
1324     //////////////////////////////////////////////////////////
1325     // Reconstruct Xi(1530)
1326     for(Int_t EN=0; EN<fEventsToMix+1; EN++){// Event buffer loop
1327
1328       for(Int_t l=0; l<(fEvt+EN)->fNTracks; l++){// Present(EN=0) and Past(EN from 1 to fEventsToMix) event track loop
1329
1330         if(EN==0) {
1331           if((fEvt+EN)->fTracks[l].fID == pTrackXi->GetID()) continue;
1332           if((fEvt+EN)->fTracks[l].fID == nTrackXi->GetID()) continue;
1333           if((fEvt+EN)->fTracks[l].fID == bTrackXi->GetID()) continue;
1334         }
1335
1336         fXiTrack->Set(xiVtx, xiP, fCovMatrix, Short_t(xiCharge));
1337
1338         if(!fESDTrack4) continue;
1339         fESDTrack4->Set((fEvt+EN)->fTracks[l].fX, (fEvt+EN)->fTracks[l].fP, (fEvt+EN)->fTracks[l].fCov, (fEvt+EN)->fTracks[l].fCharge);
1340         if(fAODcase){
1341           if((Bool_t)(((1<<5) & (fEvt+EN)->fTracks[l].fFilterMap) == 0)) continue;// AOD filterbit cut, "Standard cuts with tight dca"
1342         }else{
1343           fDecayParameters[8] = (fEvt+EN)->fTracks[l].fDCAXY;// DCA Vtx pion third
1344           if((fEvt+EN)->fTracks[l].fDCAZ > 2) continue;
1345           if( (((fEvt+EN)->fTracks[l].fStatus)&AliESDtrack::kITSrefit)==0) continue;// Require itsrefit
1346           // no Chi^2 cut applied for ESDs.  Info not available in my track structure.
1347         }
1348
1349         if(fabs((fEvt+EN)->fTracks[l].fEta) > 0.8) continue;
1350
1351         fDecayParameters[3] = (fEvt+EN)->fTracks[l].fNclusTPC;
1352
1353         AliVertex *XiStarVtx = new AliVertex((fEvt+EN)->fTracks[l].fX,0,0);
1354         //fESDTrack4->PropagateToDCA(fXiTrack, bField);// Propagate tracks to dca, both tracks are budged
1355         if(!(fXiTrack->PropagateToDCA(XiStarVtx, bField, 3))) continue;// Propagate tracks to dca, version which assumes fESDTrack4 is already primary
1356         /////////////
1357         fXiTrack->GetPxPyPz(pDaughter1);
1358         fXiTrack->GetXYZ(xDaughter1);
1359         fESDTrack4->GetPxPyPz(pDaughter2);
1360         fESDTrack4->GetXYZ(xDaughter2);
1361         //////////////////////////
1362
1363
1364
1365         //xiStarVtx[0] = (xDaughter1[0]+xDaughter2[0])/2.;
1366         //xiStarVtx[1] = (xDaughter1[1]+xDaughter2[1])/2.;
1367         //xiStarVtx[2] = (xDaughter1[2]+xDaughter2[2])/2.;
1368         //decayLength = sqrt(pow(xiStarVtx[0]-primaryVtx[0],2)+pow(xiStarVtx[1]-primaryVtx[1],2)+pow(xiStarVtx[2]-primaryVtx[2],2));
1369
1370         px1=pDaughter1[0];
1371         py1=pDaughter1[1];
1372         pz1=pDaughter1[2];
1373         px2=pDaughter2[0];
1374         py2=pDaughter2[1];
1375         pz2=pDaughter2[2];
1376
1377         p1sq=px1*px1+py1*py1+pz1*pz1;
1378         p2sq=px2*px2+py2*py2+pz2*pz2;
1379         if(p1sq <=0 || p2sq <=0) continue;
1380
1381         e1=sqrt(p1sq+fTrueMassXi*fTrueMassXi);
1382         e2=sqrt(p2sq+fTrueMassPi*fTrueMassPi);
1383         angle=px1*px2+py1*py2+pz1*pz2;
1384         xiStarMass=fTrueMassXi*fTrueMassXi+fTrueMassPi*fTrueMassPi+2.*e1*e2-2.*angle;
1385         if(xiStarMass<0.) xiStarMass=1.e-8;
1386         xiStarMass=sqrt(xiStarMass);
1387
1388
1389         xiStarP[0] = px1+px2;
1390         xiStarP[1] = py1+py2;
1391         xiStarP[2] = pz1+pz2;
1392         xiStarMom = sqrt(pow(xiStarP[0],2)+pow(xiStarP[1],2)+pow(xiStarP[2],2));
1393         if(xiStarMom==0) continue; // So one of the following lines doesnt break
1394         xiStarPt = sqrt(xiStarP[0]*xiStarP[0] + xiStarP[1]*xiStarP[1]);
1395         xiStarY = .5*log( ((e1+e2) + xiStarP[2])/((e1+e2) - xiStarP[2]));
1396         //xiStarE = e1 + e2;
1397
1398
1399         //if( (xiStarP[0]*(xiStarVtx[0]-primaryVtx[0]) + xiStarP[1]*(xiStarVtx[1]-primaryVtx[1]) + xiStarP[2]*(xiStarVtx[2]-primaryVtx[2]))/xiStarMom/decayLength < fXiStarCosTheta) continue;
1400
1401         for(int cv=0; cv<kNCutVariations; cv++){
1402           if(fDecayParameters[0] < fCutValues[cv][0]) continue;// Nclus proton
1403           if(fDecayParameters[1] < fCutValues[cv][1]) continue;// Nclus pion first
1404           if(fDecayParameters[2] < fCutValues[cv][2]) continue;// Nclus pion second
1405           if(fDecayParameters[3] < fCutValues[cv][3]) continue;// Nclus pion third
1406           //
1407           if(fDecayParameters[4] < fCutValues[cv][4]) continue;// DCAVtx proton
1408           if(fDecayParameters[5] < fCutValues[cv][5]) continue;// DCAVtx pion first
1409           if(fDecayParameters[6] < fCutValues[cv][6]) continue;// DCAVtx pion second
1410           if(fDecayParameters[7] < fCutValues[cv][7]) continue;// DCAVtx Lambda
1411           if(cv!=8) {if(fDecayParameters[8] > (0.0182 + 0.035/pow((fEvt+EN)->fTracks[l].fPt,1.01))) continue;}// DCAVtx pion third
1412           else {if(fDecayParameters[8] > fCutValues[cv][8]) continue;}// DCAVtx pion third
1413           //
1414           if(fDecayParameters[9] > fCutValues[cv][9]) continue;// DCAV proton-pion
1415           if(fDecayParameters[10] > fCutValues[cv][10]) continue;// DCAV Lambda-pion
1416           //
1417           if(fDecayParameters[11] < fCutValues[cv][11]) continue;// Rxy Lambda
1418           if(fDecayParameters[12] < fCutValues[cv][12]) continue;// Rxy Xi
1419           //
1420           if(fDecayParameters[13] < fCutValues[cv][13]) continue;// Cos PA Lambda
1421           if(fDecayParameters[14] < fCutValues[cv][14]) continue;// Cos PA Xi
1422
1423
1424           if(EN==0){
1425             if(fXiTrack->Charge() == -1 &&  fESDTrack4->Charge() == -1) CutVar[cv].fXiMinusPiMinus->Fill(xiStarPt, xiStarY, xiStarMass);
1426             else if(fXiTrack->Charge() == -1 &&  fESDTrack4->Charge() == +1) CutVar[cv].fXiMinusPiPlus->Fill(xiStarPt, xiStarY, xiStarMass);
1427             else if(fXiTrack->Charge() == +1 &&  fESDTrack4->Charge() == -1) CutVar[cv].fXiPlusPiMinus->Fill(xiStarPt, xiStarY, xiStarMass);
1428             else CutVar[cv].fXiPlusPiPlus->Fill(xiStarPt, xiStarY, xiStarMass);
1429           }else {
1430             if(fXiTrack->Charge() == -1 &&  fESDTrack4->Charge() == -1) CutVar[cv].fXiMinusPiMinusbkg->Fill(xiStarPt, xiStarY, xiStarMass);
1431             else if(fXiTrack->Charge() == -1 &&  fESDTrack4->Charge() == +1) CutVar[cv].fXiMinusPiPlusbkg->Fill(xiStarPt, xiStarY, xiStarMass);
1432             else if(fXiTrack->Charge() == +1 &&  fESDTrack4->Charge() == -1) CutVar[cv].fXiPlusPiMinusbkg->Fill(xiStarPt, xiStarY, xiStarMass);
1433             else CutVar[cv].fXiPlusPiPlusbkg->Fill(xiStarPt, xiStarY, xiStarMass);
1434           }
1435
1436
1437         /*
1438         // MC associaton AOD
1439         if(fMCcase && mcXiFilled && EN==0 && fAODcase){// AOD MC's
1440
1441           MCXiStarD2 = (AliAODMCParticle*)mcArray->At(abs((fEvt)->fTracks[l].fLabel));
1442
1443           if(abs(MCXiStarD2->GetPdgCode())==kPionCode){
1444             if(MCXi->GetMother() == MCXiStarD2->GetMother()){
1445               MCXiStar = (AliAODMCParticle*)mcArray->At(MCXi->GetMother());
1446               if(abs(MCXiStar->GetPdgCode())==kXiStarCode) {
1447
1448                 if(fXiTrack->Charge() == -1 &&  fESDTrack4->Charge() == +1) fMCrecXiMinusPiPlus->Fill(xiStarPt, xiStarY, xiStarMass);
1449                 if(fXiTrack->Charge() == +1 &&  fESDTrack4->Charge() == -1) fMCrecXiPlusPiMinus->Fill(xiStarPt, xiStarY, xiStarMass);
1450
1451               }
1452             }
1453           }
1454         }
1455         */
1456
1457         // MC associaton ESD
1458           if(fMCcase && mcXiFilled && EN==0 && !fAODcase){// ESD MC's
1459             MCXiStarD2esd = (TParticle*)mcstack->Particle(abs((fEvt)->fTracks[l].fLabel));
1460
1461             if(abs(MCXiStarD2esd->GetPdgCode())==kPionCode){
1462               if(MCXiesd->GetMother(0) == MCXiStarD2esd->GetMother(0)){
1463
1464                 MCXiStaresd = (TParticle*)mcstack->Particle(abs(MCXiesd->GetMother(0)));
1465                 if(abs(MCXiStaresd->GetPdgCode())==kXiStarCode) {
1466
1467                   if(fXiTrack->Charge() == -1 &&  fESDTrack4->Charge() == +1) CutVar[cv].fMCrecXiMinusPiPlus->Fill(xiStarPt, xiStarY, xiStarMass);
1468                   if(fXiTrack->Charge() == +1 &&  fESDTrack4->Charge() == -1) CutVar[cv].fMCrecXiPlusPiMinus->Fill(xiStarPt, xiStarY, xiStarMass);
1469
1470                 }
1471               }
1472             }
1473           }
1474
1475         }// Cut Variation loop
1476       }// 3rd pion loop
1477     }// Event mixing loop
1478
1479   }// Xi loop
1480
1481
1482
1483   // Post output data.
1484   PostData(1, fOutputList);
1485
1486 }
1487 //________________________________________________________________________
1488 void AliXiStar::Terminate(Option_t *)
1489 {
1490   cout<<"Done"<<endl;
1491 }
1492 //________________________________________________________________________
1493 Double_t AliXiStar::LinearPropagateToDCA(AliESDtrack *v, AliESDtrack *t, Double_t b) {// Adapted from AliCascadeVertexer.cxx
1494  //--------------------------------------------------------------------
1495   // This function returns the DCA between the V0 and the track
1496   //--------------------------------------------------------------------
1497
1498   Double_t alpha=t->GetAlpha(), cs1=TMath::Cos(alpha), sn1=TMath::Sin(alpha);
1499   Double_t r[3]; t->GetXYZ(r);
1500   Double_t x1=r[0], y1=r[1], z1=r[2];
1501   Double_t p[3]; t->GetPxPyPz(p);
1502   Double_t px1=p[0], py1=p[1], pz1=p[2];
1503
1504   Double_t x2[3]={0};
1505   Double_t p2[3]={0};
1506   Double_t vx2,vy2,vz2;     // position and momentum of V0
1507   Double_t px2,py2,pz2;
1508
1509   v->GetXYZ(x2);
1510   v->GetPxPyPz(p2);
1511   vx2=x2[0], vy2=x2[1], vz2=x2[2];
1512   px2=p2[0], py2=p2[1], pz2=p2[2];
1513
1514 // calculation dca
1515
1516   Double_t dd= Det(vx2-x1,vy2-y1,vz2-z1,px1,py1,pz1,px2,py2,pz2);
1517   Double_t ax= Det(py1,pz1,py2,pz2);
1518   Double_t ay=-Det(px1,pz1,px2,pz2);
1519   Double_t az= Det(px1,py1,px2,py2);
1520
1521   Double_t dca=TMath::Abs(dd)/TMath::Sqrt(ax*ax + ay*ay + az*az);
1522
1523 //points of the DCA
1524   Double_t t1 = Det(vx2-x1,vy2-y1,vz2-z1,px2,py2,pz2,ax,ay,az)/
1525                 Det(px1,py1,pz1,px2,py2,pz2,ax,ay,az);
1526
1527   x1 += px1*t1; y1 += py1*t1; //z1 += pz1*t1;
1528
1529
1530   //propagate track to the points of DCA
1531
1532   x1=x1*cs1 + y1*sn1;
1533   if (!t->PropagateTo(x1,b)) {
1534     Error("PropagateToDCA","Propagation failed !");
1535     return 1.e+33;
1536   }
1537
1538   return dca;
1539 }
1540
1541
1542 //________________________________________________________________________
1543 Double_t AliXiStar::Det(Double_t a00, Double_t a01, Double_t a10, Double_t a11) const {// Taken from AliCascadeVertexer
1544   //--------------------------------------------------------------------
1545   // This function calculates locally a 2x2 determinant
1546   //--------------------------------------------------------------------
1547   return a00*a11 - a01*a10;
1548 }
1549 //________________________________________________________________________
1550 Double_t AliXiStar::Det(Double_t a00,Double_t a01,Double_t a02,
1551                                  Double_t a10,Double_t a11,Double_t a12,
1552                               Double_t a20,Double_t a21,Double_t a22) const {// Taken from AliCascadeVertexer
1553   //--------------------------------------------------------------------
1554   // This function calculates locally a 3x3 determinant
1555   //--------------------------------------------------------------------
1556   return  a00*Det(a11,a12,a21,a22)-a01*Det(a10,a12,a20,a22)+a02*Det(a10,a11,a20,a21);
1557 }
1558
1559
1560