PWGLF/STRANGENESS/LambdaK0PbPb/AliAnalysisTaskCTauPbPb.cxx

   1 #include <TCanvas.h>
   2 #include <TTree.h>
   3 #include <TFile.h>
   4 #include <TH1F.h>
   5 #include <TH2F.h>
   6 #include <TH3F.h>
   7 #include <TPDGCode.h>
   8 #include <TDatabasePDG.h>
   9 #include <TParticlePDG.h>
  10 #include <TParticle.h>
  11 #include <TROOT.h>
  12
  13 #include "AliESDEvent.h"
  14 #include "AliESDv0.h"
  15 #include "AliESDcascade.h"
  16
  17 #include "AliCentrality.h"
  18
  19 #include "AliMCEvent.h"
  20 #include "AliStack.h"
  21
  22 #include "AliPID.h"
  23 #include "AliPIDResponse.h"
  24
  25 #include "AliInputEventHandler.h"
  26 #include "AliAnalysisManager.h"
  27
  28 #include "AliAnalysisTaskCTauPbPb.h"
  29
  30 extern TROOT *gROOT;
  31
  32 ClassImp(AliAnalysisTaskCTauPbPb)
  33
  34
  35 //
  36 //  This is a little task for checking the c*tau of the strange particles
  37 //
  38
  39 AliAnalysisTaskCTauPbPb::AliAnalysisTaskCTauPbPb(const char *name) :
  40 AliAnalysisTaskSE(name),
  41 fIsMC(kFALSE),
  42 fCMin(0.),
  43 fCMax(90.),
  44 fOutput(0),
  45 fMult(0),
  46 fdEdx(0),
  47 fdEdxPid(0),
  48
  49 fK0sM(0),
  50 fK0sSi(0),
  51 fK0sMC(0),
  52 fK0sAs(0),
  53
  54 fLambdaM(0),
  55 fLambdaSi(0),
  56 fLambdaMC(0),
  57 fLambdaAs(0),
  58
  59 fLambdaEff(0),
  60 fLambdaPt(0),
  61
  62 fLambdaBarM(0),
  63 fLambdaBarSi(0),
  64 fLambdaBarMC(0),
  65 fLambdaBarAs(0),
  66
  67 fLambdaBarEff(0),
  68 fLambdaBarPt(0),
  69
  70 fLambdaFromXi(0),
  71 fXiM(0),
  72 fXiSiP(0),
  73
  74 fLambdaBarFromXiBar(0),
  75 fXiBarM(0),
  76 fXiBarSiP(0)
  77 {
  78   // Constructor. Initialization of pointers
  79   DefineOutput(1, TList::Class());
  80 }
  81
  82 void AliAnalysisTaskCTauPbPb::UserCreateOutputObjects()
  83 {
  84   Int_t    nbins=100;  // number of bins
  85   Double_t ltMax=100.;
  86   Double_t ptMax=10.;
  87
  88   fOutput = new TList();
  89   fOutput->SetOwner();
  90
  91
  92   fMult=new TH1F("fMult","Multiplicity",1100,0.,3300);
  93   fMult->GetXaxis()->SetTitle("N tracks");
  94   fOutput->Add(fMult);
  95
  96   fdEdx=new TH2F("fdEdx","dE/dx",50,0.2,3,50,0.,6.);
  97   fOutput->Add(fdEdx);
  98
  99   fdEdxPid=new TH2F("fdEdxPid","dE/dx with PID",50,0.2,3,50,0.,6.);
 100   fOutput->Add(fdEdxPid);
 101
 102   fK0sM =
 103   new TH2F("fK0sM", "Mass for K^{0}_{s}", nbins/2, 0.448, 0.548, nbins,0.,ptMax);
 104   fK0sM->GetXaxis()->SetTitle("Mass [GeV/c]");
 105   fOutput->Add(fK0sM);
 106
 107   fK0sSi =
 108   new TH2F("fK0sSi","L_{T} vs p_{T} for K^{0}_{s}, side-band subtracted",
 109   nbins,0.,ptMax,nbins,0.,ltMax);
 110   fK0sSi->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 111   fK0sSi->GetYaxis()->SetTitle("L_{T} [cm]");
 112   fOutput->Add(fK0sSi);
 113
 114   fK0sMC =
 115   new TH2F("fK0sMC","L_{T} vs p_{T} for K^{0}_{s}, from MC stack",
 116   nbins,0.,ptMax,nbins,0.,ltMax);
 117   fK0sMC->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 118   fK0sMC->GetYaxis()->SetTitle("L_{T} [cm]");
 119   fOutput->Add(fK0sMC);
 120
 121   fK0sAs =
 122   new TH2F("fK0sAs", "L_{T} vs p_{T} for K^{0}_{s}, associated",
 123   nbins,0.,ptMax,nbins,0.,ltMax);
 124   fK0sAs->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 125   fK0sAs->GetYaxis()->SetTitle("L_{T} [cm]");
 126   fOutput->Add(fK0sAs);
 127
 128   //----------------------
 129
 130   fLambdaM =
 131   new TH2F("fLambdaM","Mass for \\Lambda", nbins, 1.065, 1.165,nbins,0.,ptMax);
 132   fLambdaM->GetXaxis()->SetTitle("Mass [GeV/c]");
 133   fOutput->Add(fLambdaM);
 134
 135   fLambdaSi =
 136   new TH2F("fLambdaSi","L_{T} vs p_{T} for \\Lambda, side-band subtructed",
 137   nbins,0.,ptMax,nbins,0.,ltMax);
 138   fLambdaSi->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 139   fLambdaSi->GetYaxis()->SetTitle("L_{T} [cm]");
 140   fOutput->Add(fLambdaSi);
 141
 142   fLambdaMC =
 143   new TH2F("fLambdaMC","c\\tau for \\Lambda, from MC stack",
 144   nbins,0.,ptMax,nbins,0.,ltMax);
 145   fLambdaMC->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 146   fLambdaMC->GetYaxis()->SetTitle("L_{T} [cm]");
 147   fOutput->Add(fLambdaMC);
 148
 149   fLambdaAs =
 150   new TH2F("fLambdaAs","c\\tau for \\Lambda, associated",
 151   nbins,0.,ptMax,nbins,0.,ltMax);
 152   fLambdaAs->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 153   fLambdaAs->GetYaxis()->SetTitle("L_{T} [cm]");
 154   fOutput->Add(fLambdaAs);
 155
 156   //----------------------
 157
 158   fLambdaEff=fLambdaAs->ProjectionX();
 159   fLambdaEff->SetName("fLambdaEff");
 160   fLambdaEff->SetTitle("Efficiency for #Lambda");
 161   fOutput->Add(fLambdaEff);
 162
 163   fLambdaPt=fLambdaAs->ProjectionX();
 164   fLambdaPt->SetName("fLambdaPt");
 165   fLambdaPt->SetTitle("Raw #Lambda pT spectrum");
 166   fOutput->Add(fLambdaPt);
 167
 168   //----------------------
 169
 170   fLambdaBarM =
 171   new TH2F("fLambdaBarM","Mass for anti-\\Lambda", nbins, 1.065, 1.165,nbins,0.,ptMax);
 172   fLambdaBarM->GetXaxis()->SetTitle("Mass [GeV/c]");
 173   fOutput->Add(fLambdaBarM);
 174
 175   fLambdaBarSi =
 176   new TH2F("fLambdaBarSi","L_{T} vs p_{T} for anti-\\Lambda, side-band subtructed",
 177   nbins,0.,ptMax,nbins,0.,ltMax);
 178   fLambdaBarSi->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 179   fLambdaBarSi->GetYaxis()->SetTitle("L_{T} [cm]");
 180   fOutput->Add(fLambdaBarSi);
 181
 182   fLambdaBarMC =
 183   new TH2F("fLambdaBarMC","c\\tau for anti-\\Lambda, from MC stack",
 184   nbins,0.,ptMax,nbins,0.,ltMax);
 185   fLambdaBarMC->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 186   fLambdaBarMC->GetYaxis()->SetTitle("L_{T} [cm]");
 187   fOutput->Add(fLambdaBarMC);
 188
 189   fLambdaBarAs =
 190   new TH2F("fLambdaBarAs","c\\tau for anti-\\Lambda, associated",
 191   nbins,0.,ptMax,nbins,0.,ltMax);
 192   fLambdaBarAs->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 193   fLambdaBarAs->GetYaxis()->SetTitle("L_{T} [cm]");
 194   fOutput->Add(fLambdaBarAs);
 195
 196
 197   //----------------------
 198
 199   fLambdaBarEff=fLambdaBarAs->ProjectionX();
 200   fLambdaBarEff->SetName("fLambdaBarEff");
 201   fLambdaBarEff->SetTitle("Efficiency for anti-#Lambda");
 202   fOutput->Add(fLambdaBarEff);
 203
 204   fLambdaBarPt=fLambdaBarAs->ProjectionX();
 205   fLambdaBarPt->SetName("fLambdaBarPt");
 206   fLambdaBarPt->SetTitle("Raw anti-#Lambda pT spectrum");
 207   fOutput->Add(fLambdaBarPt);
 208
 209   //----------------------
 210
 211   fLambdaFromXi=new TH3F("fLambdaFromXi","L_{T} vs p_{T} vs p_{T} of \\Xi for \\Lambda from \\Xi",
 212   nbins,0.,ptMax,nbins,0.,ltMax,33,0.,ptMax+2);
 213   fOutput->Add(fLambdaFromXi);
 214
 215   fXiM  =
 216   new TH2F("fXiM", "\\Xi mass distribution", 50, 1.271, 1.371,33,0.,ptMax+2);
 217   fOutput->Add(fXiM);
 218
 219   fXiSiP  = new TH1F("fXiSiP", "Pt for \\Xi, side-band subracted",
 220   33,0.,ptMax+2);
 221   fOutput->Add(fXiSiP);
 222
 223
 224   fLambdaBarFromXiBar=new TH3F("fLambdaBarFromXiBar","L_{T} vs p_{T} vs p_{T} of anti-\\Xi for anti-\\Lambda from anti-\\Xi",
 225   nbins,0.,ptMax,nbins,0.,ltMax,33,0.,ptMax+2);
 226   fOutput->Add(fLambdaBarFromXiBar);
 227
 228   fXiBarM  =
 229   new TH2F("fXiBarM", "anti-\\Xi mass distribution", 50, 1.271, 1.371,33,0.,ptMax+2);
 230   fOutput->Add(fXiBarM);
 231
 232   fXiBarSiP  = new TH1F("fXiBarSiP", "Pt for anti-\\Xi, side-band subracted",
 233   33,0.,ptMax+2);
 234   fOutput->Add(fXiBarSiP);
 235
 236
 237   PostData(1, fOutput);
 238 }
 239
 240 static Bool_t AcceptTrack(const AliESDtrack *t) {
 241   if (!t->IsOn(AliESDtrack::kTPCrefit)) return kFALSE;
 242   if (t->GetKinkIndex(0)>0) return kFALSE;
 243
 244   Float_t nCrossedRowsTPC = t->GetTPCClusterInfo(2,1);
 245   if (nCrossedRowsTPC < 70) return kFALSE;
 246   Int_t findable=t->GetTPCNclsF();
 247   if (findable <= 0) return kFALSE;
 248   if (nCrossedRowsTPC/findable < 0.8) return kFALSE;
 249
 250   if (TMath::Abs(t->Eta()) > 0.8) return kFALSE;
 251
 252   return kTRUE;
 253 }
 254
 255 static Bool_t AcceptV0(const AliESDv0 *v0, const AliESDEvent *esd) {
 256
 257   if (v0->GetOnFlyStatus()) return kFALSE;
 258
 259   Int_t nidx=TMath::Abs(v0->GetNindex());
 260   AliESDtrack *ntrack=esd->GetTrack(nidx);
 261   if (!AcceptTrack(ntrack)) return kFALSE;
 262
 263   Int_t pidx=TMath::Abs(v0->GetPindex());
 264   AliESDtrack *ptrack=esd->GetTrack(pidx);
 265   if (!AcceptTrack(ptrack)) return kFALSE;
 266
 267   Float_t xy,z0;
 268   ntrack->GetImpactParameters(xy,z0);
 269   if (TMath::Abs(xy)<0.1) return kFALSE;
 270   ptrack->GetImpactParameters(xy,z0);
 271   if (TMath::Abs(xy)<0.1) return kFALSE;
 272
 273   Double_t dca=v0->GetDcaV0Daughters();
 274   if (dca>1.0) return kFALSE;
 275   //if (dca>0.7) return kFALSE;
 276   //if (dca>0.4) return kFALSE;
 277
 278   Double_t cpa=v0->GetV0CosineOfPointingAngle();
 279   if (cpa<0.998) return kFALSE;
 280   //if (cpa<0.99875) return kFALSE;
 281   //if (cpa<0.9995) return kFALSE;
 282
 283   Double_t xx,yy,zz; v0->GetXYZ(xx,yy,zz);
 284   Double_t r2=xx*xx + yy*yy;
 285   if (r2<0.9*0.9) return kFALSE;
 286   if (r2>100*100) return kFALSE;
 287
 288   return kTRUE;
 289 }
 290
 291 static Bool_t AcceptCascade(const AliESDcascade *cs, const AliESDEvent *esd) {
 292
 293   Int_t bidx=TMath::Abs(cs->GetBindex());
 294   AliESDtrack *btrack=esd->GetTrack(bidx);
 295   if (!AcceptTrack(btrack)) return kFALSE;
 296
 297   Float_t xy,z0;
 298   btrack->GetImpactParameters(xy,z0);
 299   if (TMath::Abs(xy)<0.03) return kFALSE;
 300
 301   const AliESDVertex *vtx=esd->GetPrimaryVertexSPD();
 302   if (!vtx->GetStatus()) {
 303      vtx=esd->GetPrimaryVertexTracks();
 304      if (!vtx->GetStatus()) return kFALSE;
 305   }
 306   Double_t xv=vtx->GetXv(), yv=vtx->GetYv(), zv=vtx->GetZv();
 307   if (cs->GetCascadeCosineOfPointingAngle(xv,yv,zv) < 0.999) return kFALSE;
 308
 309   if (cs->GetDcaXiDaughters() > 0.3) return kFALSE;
 310
 311   return kTRUE;
 312 }
 313
 314 static Bool_t AcceptPID(const AliPIDResponse *pidResponse,
 315                         const AliESDtrack *ptrack, AliStack *stack) {
 316
 317   const AliExternalTrackParam *par=ptrack->GetInnerParam();
 318   if (!par) return kTRUE;
 319   if (par->GetP() > 1.) return kTRUE;
 320
 321
 322   if (stack) {
 323     // MC PID
 324     Int_t plab=TMath::Abs(ptrack->GetLabel());
 325     TParticle *pp=stack->Particle(plab);
 326     if (!pp) return kTRUE;
 327     if (pp->GetPDG()->Charge() > 0) {
 328        if (pp->GetPdgCode() == kProton)    return kTRUE;
 329     } else {
 330        if (pp->GetPdgCode() == kProtonBar) return kTRUE;
 331     }
 332   } else {
 333     // Real PID
 334     Double_t nsig=pidResponse->NumberOfSigmasTPC(ptrack,AliPID::kProton);
 335     if (TMath::Abs(nsig) < 3.) return kTRUE;
 336   }
 337
 338   return kFALSE;
 339 }
 340
 341 static TParticle*
 342 AssociateV0(const AliESDtrack *ptrack,const AliESDtrack *ntrack,AliStack *stack,
 343 TParticle *&mcp) {
 344   //
 345   // Try to associate the V0 with the daughters ptrack and ntrack
 346   // with the Monte Carlo
 347   //
 348   if (!stack) return 0;
 349
 350   Int_t nlab=TMath::Abs(ntrack->GetLabel());
 351   TParticle *n=stack->Particle(nlab);
 352   if (!n) return 0;
 353
 354   Int_t plab=TMath::Abs(ptrack->GetLabel());
 355   TParticle *p=stack->Particle(plab);
 356   if (!p) return 0;
 357
 358   Int_t imp=p->GetFirstMother();
 359   if (imp<0) return 0;
 360   if (imp>=stack->GetNtrack()) return 0;
 361   TParticle *p0=stack->Particle(imp); // V0 particle, mother of the pos. track
 362   if (!p0) return 0;
 363
 364   Int_t imn=n->GetFirstMother();
 365   if (imp != imn) {  // Check decays of the daughters
 366     return 0; // Fixme
 367   }
 368
 369   Int_t code=p0->GetPdgCode();
 370   if (code != kK0Short)
 371      if (code != kLambda0)
 372         if (code != kLambda0Bar) return 0;
 373
 374   mcp=p;
 375
 376   return p0;
 377 }
 378
 379
 380 void AliAnalysisTaskCTauPbPb::UserExec(Option_t *)
 381 {
 382   const Double_t yMax=0.5;
 383   const Double_t pMin=0.0;
 384   const Double_t lMax=0.001;
 385
 386   AliESDEvent *esd=(AliESDEvent *)InputEvent();
 387
 388   if (!esd) {
 389     Printf("ERROR: esd not available");
 390     return;
 391   }
 392
 393   // Vertex selection
 394   const AliESDVertex *vtx=esd->GetPrimaryVertexSPD();
 395   if (!vtx->GetStatus()) {
 396      vtx=esd->GetPrimaryVertexTracks();
 397      if (!vtx->GetStatus()) return;
 398   }
 399   Double_t xv=vtx->GetXv(), yv=vtx->GetYv(), zv=vtx->GetZv();
 400
 401   if (TMath::Abs(zv) > 10.) return ;
 402
 403
 404   // Physics selection
 405   AliAnalysisManager *mgr= AliAnalysisManager::GetAnalysisManager();
 406   AliInputEventHandler *hdr=(AliInputEventHandler*)mgr->GetInputEventHandler();
 407   UInt_t maskIsSelected = hdr->IsEventSelected();
 408   Bool_t isSelected = (maskIsSelected & AliVEvent::kMB);
 409   if (!isSelected) return;
 410
 411
 412   fMult->Fill(-100); //event counter
 413
 414
 415   // Centrality selection
 416   AliCentrality *cent=esd->GetCentrality();
 417   if (!cent->IsEventInCentralityClass(fCMin,fCMax,"V0M")) return;
 418
 419   AliPIDResponse *pidResponse = hdr->GetPIDResponse();
 420
 421
 422   //+++++++ MC
 423   AliStack *stack = 0x0;
 424   Double_t mcXv=0., mcYv=0., mcZv=0.;
 425
 426   if (fIsMC) {
 427      AliMCEvent *mcEvent = MCEvent();
 428      stack = mcEvent->Stack();
 429      if (!stack) {
 430         Printf("ERROR: stack not available");
 431         return;
 432      }
 433
 434      const AliVVertex *mcVtx=mcEvent->GetPrimaryVertex();
 435
 436      mcXv=mcVtx->GetX(); mcYv=mcVtx->GetY(); mcZv=mcVtx->GetZ();
 437
 438      Int_t ntrk=stack->GetNtrack();
 439      while (ntrk--) {
 440        TParticle *p0=stack->Particle(ntrk);
 441        Int_t code=p0->GetPdgCode();
 442        if (code != kK0Short)
 443          if (code != kLambda0)
 444             if (code != kLambda0Bar) continue;
 445
 446        Int_t plab=p0->GetFirstDaughter(), nlab=p0->GetLastDaughter();
 447        if (nlab==plab) continue;
 448        TParticle *part = stack->Particle(plab);
 449        if (!part) continue;
 450        TParticlePDG *partPDG = part->GetPDG();
 451        if (!partPDG) continue;
 452        Double_t charge=partPDG->Charge();
 453        if (charge == 0.) continue;
 454
 455        Double_t pt=p0->Pt();
 456        if (pt<pMin) continue;
 457        if (TMath::Abs(p0->Y())>yMax) continue;
 458
 459        Double_t x=p0->Vx(), y=p0->Vy(), z=p0->Vz();
 460        Double_t dx=mcXv-x, dy=mcYv-y, dz=mcZv-z;
 461        Double_t l=TMath::Sqrt(dx*dx + dy*dy + dz*dz);
 462
 463        if (l > lMax) continue; // secondary V0
 464
 465        x=part->Vx(); y=part->Vy();
 466        dx=mcXv-x; dy=mcYv-y;
 467        Double_t lt=TMath::Sqrt(dx*dx + dy*dy);
 468
 469        switch (code) {
 470        case kK0Short:
 471           fK0sMC->Fill(pt,lt);
 472           break;
 473        case kLambda0:
 474           fLambdaMC->Fill(pt,lt);
 475           break;
 476        case kLambda0Bar:
 477           fLambdaBarMC->Fill(pt,lt);
 478           break;
 479        default: break;
 480        }
 481      }
 482   }
 483   //+++++++
 484
 485
 486   Int_t ntrk1=esd->GetNumberOfTracks();
 487   Int_t mult=0;
 488   for (Int_t i=0; i<ntrk1; i++) {
 489     AliESDtrack *t=esd->GetTrack(i);
 490     if (!t->IsOn(AliESDtrack::kTPCrefit)) continue;
 491     Float_t xy,z0;
 492     t->GetImpactParameters(xy,z0);
 493     if (TMath::Abs(xy)>3.) continue;
 494     if (TMath::Abs(z0)>3.) continue;
 495     Double_t pt=t->Pt(),pz=t->Pz();
 496     if (TMath::Abs(pz/pt)>0.8) continue;
 497     mult++;
 498
 499     Double_t p=t->GetInnerParam()->GetP();
 500     Double_t dedx=t->GetTPCsignal()/47.;
 501     fdEdx->Fill(p,dedx,1);
 502
 503     Double_t nsig=pidResponse->NumberOfSigmasTPC(t,AliPID::kProton);
 504     if (TMath::Abs(nsig) < 3.) fdEdxPid->Fill(p,dedx,1);
 505
 506   }
 507   fMult->Fill(mult);
 508
 509
 510   Int_t nv0 = esd->GetNumberOfV0s();
 511   while (nv0--) {
 512       AliESDv0 *v0=esd->GetV0(nv0);
 513
 514       Double_t pt=v0->Pt();
 515       if (pt < pMin) continue;
 516
 517       if (!AcceptV0(v0,esd)) continue;
 518
 519       Int_t nidx=TMath::Abs(v0->GetNindex());
 520       AliESDtrack *ntrack=esd->GetTrack(nidx);
 521       Int_t pidx=TMath::Abs(v0->GetPindex());
 522       AliESDtrack *ptrack=esd->GetTrack(pidx);
 523
 524       Double_t x,y,z; v0->GetXYZ(x,y,z);
 525       Double_t dx1=x-xv, dy1=y-yv;
 526       Double_t lt=TMath::Sqrt(dx1*dx1 + dy1*dy1);
 527
 528       if (lt/pt > 3*7.89/1.1157) continue;
 529
 530       //--- V0 switches
 531       Bool_t isK0s=kTRUE;
 532       Bool_t isLambda=kTRUE;
 533       Bool_t isLambdaBar=kTRUE;
 534
 535       if (0.4977*lt/pt > 3*2.68) isK0s=kFALSE;
 536       if (1.1157*lt/pt > 3*7.89) isLambdaBar=isLambda=kFALSE;
 537
 538       if (!AcceptPID(pidResponse, ptrack, stack)) isLambda=kFALSE;
 539       if (!AcceptPID(pidResponse, ntrack, stack)) isLambdaBar=kFALSE;
 540
 541       Double_t yK0s=TMath::Abs(v0->RapK0Short());
 542       Double_t yLam=TMath::Abs(v0->RapLambda());
 543       if (yK0s > yMax) isK0s=kFALSE;
 544       if (yLam > yMax) isLambda=isLambdaBar=kFALSE;
 545       //---
 546
 547       Double_t mass=0., m=0., s=0.;
 548       if (isK0s) {
 549          v0->ChangeMassHypothesis(kK0Short);
 550
 551          mass=v0->GetEffMass();
 552          fK0sM->Fill(mass,pt);
 553
 554          m=TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass();
 555          s=0.0044 + (0.008-0.0044)/(10-1)*(pt - 1.);
 556          if (TMath::Abs(m-mass) < 3*s) {
 557             fK0sSi->Fill(pt,lt);
 558          } else {
 559             isK0s=kFALSE;
 560          }
 561          if (TMath::Abs(m-mass + 4.5*s) < 1.5*s) {
 562             fK0sSi->Fill(pt,lt,-1);
 563          }
 564          if (TMath::Abs(m-mass - 4.5*s) < 1.5*s) {
 565             fK0sSi->Fill(pt,lt,-1);
 566          }
 567       }
 568
 569       if (isLambda) {
 570          v0->ChangeMassHypothesis(kLambda0);
 571
 572          mass=v0->GetEffMass();
 573          fLambdaM->Fill(mass,pt);
 574
 575          m=TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass();
 576          //s=0.0027 + (0.004-0.0027)/(10-1)*(pt-1);
 577          //s=0.0015 + (0.002-0.0015)/(2.6-1)*(pt-1);
 578          s=0.0023 + (0.004-0.0023)/(6-1)*(pt-1);
 579          if (TMath::Abs(m-mass) < 3*s) {
 580             fLambdaSi->Fill(pt,lt);
 581          } else {
 582             isLambda=kFALSE;
 583          }
 584          if (TMath::Abs(m-mass + 4.5*s) < 1.5*s) {
 585             fLambdaSi->Fill(pt,lt,-1);
 586          }
 587          if (TMath::Abs(m-mass - 4.5*s) < 1.5*s) {
 588             fLambdaSi->Fill(pt,lt,-1);
 589          }
 590       }
 591
 592       if (isLambdaBar) {
 593          v0->ChangeMassHypothesis(kLambda0Bar);
 594
 595          mass=v0->GetEffMass();
 596          fLambdaBarM->Fill(mass,pt);
 597
 598          m=TDatabasePDG::Instance()->GetParticle(kLambda0Bar)->Mass();
 599          //s=0.0027 + (0.004-0.0027)/(10-1)*(pt-1);
 600          //s=0.0015 + (0.002-0.0015)/(2.6-1)*(pt-1);
 601          s=0.0023 + (0.004-0.0023)/(6-1)*(pt-1);
 602          if (TMath::Abs(m-mass) < 3*s) {
 603             fLambdaBarSi->Fill(pt,lt);
 604          } else {
 605             isLambdaBar=kFALSE;
 606          }
 607          if (TMath::Abs(m-mass + 4.5*s) < 1.5*s) {
 608             fLambdaBarSi->Fill(pt,lt,-1);
 609          }
 610          if (TMath::Abs(m-mass - 4.5*s) < 1.5*s) {
 611             fLambdaBarSi->Fill(pt,lt,-1);
 612          }
 613       }
 614
 615       if (!fIsMC) continue;
 616
 617       //++++++ MC
 618       if (!isK0s)
 619          if (!isLambda)
 620              if (!isLambdaBar) continue;//check MC only for the accepted V0s
 621
 622       TParticle *mcp=0;
 623       TParticle *mc0=AssociateV0(ptrack,ntrack,stack,mcp);
 624       if (!mc0) continue;
 625
 626       Double_t ptAs=mc0->Pt();
 627       if (ptAs < pMin) continue;
 628       Double_t yAs=mc0->Y();
 629       if (TMath::Abs(yAs) > yMax) continue;
 630
 631       Int_t code=mc0->GetPdgCode();
 632
 633       Double_t dx = mcXv - mcp->Vx(), dy = mcYv - mcp->Vy();
 634       Double_t ltAs=TMath::Sqrt(dx*dx + dy*dy);
 635
 636       Double_t dz=mcZv - mc0->Vz(); dy=mcYv - mc0->Vy(); dx=mcXv - mc0->Vx();
 637       Double_t l = TMath::Sqrt(dx*dx + dy*dy + dz*dz);
 638       if (l<lMax) { // Primary V0s
 639          switch (code) {
 640          case kK0Short:
 641             if (isK0s)       fK0sAs->Fill(ptAs,ltAs);
 642             break;
 643          case kLambda0:
 644             if (isLambda)    fLambdaAs->Fill(ptAs,ltAs);
 645             break;
 646          case kLambda0Bar:
 647             if (isLambdaBar) fLambdaBarAs->Fill(ptAs,ltAs);
 648             break;
 649          default: break;
 650          }
 651       } else {
 652          if (code==kK0Short) continue;
 653
 654          Int_t nx=mc0->GetFirstMother();
 655          if (nx<0) continue;
 656          if (nx>=stack->GetNtrack()) continue;
 657          TParticle *xi=stack->Particle(nx);
 658          if (!xi) continue;
 659          Int_t xcode=xi->GetPdgCode();
 660
 661          switch (code) {
 662          case kLambda0:
 663             if (isLambda)
 664             if ((xcode==kXiMinus) || (xcode==3322))
 665                fLambdaFromXi->Fill(ptAs,ltAs,xi->Pt());
 666             break;
 667          case kLambda0Bar:
 668             if (isLambdaBar)
 669             if ((xcode==kXiPlusBar)||(xcode==-3322))
 670                fLambdaBarFromXiBar->Fill(ptAs,ltAs,xi->Pt());
 671             break;
 672          default: break;
 673          }
 674       }
 675       //++++++
 676
 677   }
 678
 679   Double_t kine0;
 680   Int_t ncs=esd->GetNumberOfCascades();
 681   for (Int_t i=0; i<ncs; i++) {
 682       AliESDcascade *cs=esd->GetCascade(i);
 683
 684       Double_t pt=cs->Pt();
 685       if (pt < pMin) continue;
 686       if (TMath::Abs(cs->RapXi()) > yMax) continue;
 687       if (!AcceptCascade(cs,esd)) continue;
 688
 689       AliESDv0 *v0 = (AliESDv0*)cs;
 690       if (v0->Pt() < pMin) continue;
 691       if (TMath::Abs(v0->RapLambda()) > yMax) continue;
 692       if (!AcceptV0(v0,esd)) continue;
 693
 694       //--- Cascade switches
 695       Bool_t isXiMinus=kTRUE;
 696       Bool_t isXiPlusBar=kTRUE;
 697
 698       Int_t pidx=TMath::Abs(v0->GetPindex());
 699       AliESDtrack *ptrack=esd->GetTrack(pidx);
 700       if (!AcceptPID(pidResponse, ptrack, stack)) isXiMinus=kFALSE;
 701
 702       Int_t nidx=TMath::Abs(v0->GetNindex());
 703       AliESDtrack *ntrack=esd->GetTrack(nidx);
 704       if (!AcceptPID(pidResponse, ntrack, stack)) isXiPlusBar=kFALSE;
 705
 706       Int_t charge=cs->Charge();
 707       if (charge > 0) isXiMinus=kFALSE;
 708       if (charge < 0) isXiPlusBar=kFALSE;
 709       //---
 710
 711       if (isXiMinus) {
 712          cs->ChangeMassHypothesis(kine0,kXiMinus);
 713          Double_t mass=cs->GetEffMassXi();
 714          fXiM->Fill(mass,pt);
 715          Double_t m=TDatabasePDG::Instance()->GetParticle(kXiMinus)->Mass();
 716          //Double_t s=0.0037;
 717          Double_t s=0.002 + (0.0032-0.002)/(6-1.5)*(pt-1.5);
 718          if (TMath::Abs(m-mass) < 3*s) {
 719             fXiSiP->Fill(pt);
 720          } else {
 721             isXiMinus=kFALSE;
 722          }
 723          if (TMath::Abs(m-mass + 4.5*s) < 1.5*s) {
 724             fXiSiP->Fill(pt,-1);
 725          }
 726          if (TMath::Abs(m-mass - 4.5*s) < 1.5*s) {
 727             fXiSiP->Fill(pt,-1);
 728          }
 729       }
 730
 731       if (isXiPlusBar) {
 732          cs->ChangeMassHypothesis(kine0,kXiPlusBar);
 733          Double_t mass=cs->GetEffMassXi();
 734          fXiBarM->Fill(mass,pt);
 735          Double_t m=TDatabasePDG::Instance()->GetParticle(kXiPlusBar)->Mass();
 736          //Double_t s=0.0037;
 737          Double_t s=0.002 + (0.0032-0.002)/(6-1.5)*(pt-1.5);
 738          if (TMath::Abs(m-mass) < 3*s) {
 739             fXiBarSiP->Fill(pt);
 740          } else {
 741             isXiPlusBar=kFALSE;
 742          }
 743          if (TMath::Abs(m-mass + 4.5*s) < 1.5*s) {
 744             fXiBarSiP->Fill(pt,-1);
 745          }
 746          if (TMath::Abs(m-mass - 4.5*s) < 1.5*s) {
 747             fXiBarSiP->Fill(pt,-1);
 748          }
 749       }
 750
 751       if (!fIsMC) continue;
 752
 753       //++++++ MC
 754       if (!isXiMinus)
 755          if (!isXiPlusBar) continue;//check MC only for the accepted cascades
 756       // Here is the future association with MC
 757   }
 758
 759 }
 760
 761 void AliAnalysisTaskCTauPbPb::Terminate(Option_t *)
 762 {
 763    // The Terminate() function is the last function to be called during
 764    // a query. It always runs on the client, it can be used to present
 765    // the results graphically or save the results to file.
 766
 767   fOutput=(TList*)GetOutputData(1);
 768   if (!fOutput) {
 769      Printf("ERROR: fOutput not available");
 770      return;
 771   }
 772
 773   fMult = dynamic_cast<TH1F*>(fOutput->FindObject("fMult")) ;
 774   if (!fMult) {
 775      Printf("ERROR: fMult not available");
 776      return;
 777   }
 778
 779   fdEdx = dynamic_cast<TH2F*>(fOutput->FindObject("fdEdx")) ;
 780   if (!fdEdx) {
 781      Printf("ERROR: fdEdx not available");
 782      return;
 783   }
 784
 785   fdEdxPid = dynamic_cast<TH2F*>(fOutput->FindObject("fdEdxPid")) ;
 786   if (!fdEdxPid) {
 787      Printf("ERROR: fdEdxPid not available");
 788      return;
 789   }
 790
 791
 792   fK0sMC = dynamic_cast<TH2F*>(fOutput->FindObject("fK0sMC")) ;
 793   if (!fK0sMC) {
 794      Printf("ERROR: fK0sMC not available");
 795      return;
 796   }
 797   TH1D *k0sMcPx=fK0sMC->ProjectionX(); k0sMcPx->Sumw2();
 798   fK0sAs = dynamic_cast<TH2F*>(fOutput->FindObject("fK0sAs")) ;
 799   if (!fK0sAs) {
 800      Printf("ERROR: fK0sAs not available");
 801      return;
 802   }
 803   TH1D *k0sAsPx=fK0sAs->ProjectionX();
 804   k0sAsPx->Sumw2(); //k0sAsPx->Scale(0.69);
 805
 806
 807
 808   // Lambda histograms
 809   fLambdaFromXi = dynamic_cast<TH3F*>(fOutput->FindObject("fLambdaFromXi")) ;
 810   if (!fLambdaFromXi) {
 811      Printf("ERROR: fLambdaFromXi not available");
 812      return;
 813   }
 814   TH1D *lambdaFromXiPx=fLambdaFromXi->ProjectionX(); lambdaFromXiPx->Sumw2();
 815
 816   fLambdaMC = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaMC")) ;
 817   if (!fLambdaMC) {
 818      Printf("ERROR: fLambdaMC not available");
 819      return;
 820   }
 821   TH1D *lambdaMcPx=fLambdaMC->ProjectionX(); lambdaMcPx->Sumw2();
 822
 823   fLambdaAs = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaAs")) ;
 824   if (!fLambdaAs) {
 825      Printf("ERROR: fLambdaAs not available");
 826      return;
 827   }
 828   TH1D *lambdaAsPx=fLambdaAs->ProjectionX();
 829   lambdaAsPx->Sumw2(); //lambdaAsPx->Scale(0.64);
 830
 831   fLambdaSi = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaSi")) ;
 832   if (!fLambdaSi) {
 833      Printf("ERROR: fLambdaSi not available");
 834      return;
 835   }
 836   TH1D *lambdaSiPx=fLambdaSi->ProjectionX();
 837   lambdaSiPx->SetName("fLambdaPt");
 838   lambdaSiPx->Sumw2();
 839
 840   fLambdaEff = dynamic_cast<TH1D*>(fOutput->FindObject("fLambdaEff")) ;
 841   if (!fLambdaEff) {
 842      Printf("ERROR: fLambdaEff not available");
 843      return;
 844   }
 845   fLambdaPt = dynamic_cast<TH1D*>(fOutput->FindObject("fLambdaPt")) ;
 846   if (!fLambdaPt) {
 847      Printf("ERROR: fLambdaPt not available");
 848      return;
 849   }
 850
 851
 852   // anti-Lambda histograms
 853   fLambdaBarFromXiBar =
 854   dynamic_cast<TH3F*>(fOutput->FindObject("fLambdaBarFromXiBar")) ;
 855   if (!fLambdaBarFromXiBar) {
 856      Printf("ERROR: fLambdaBarFromXiBar not available");
 857      return;
 858   }
 859   TH1D *lambdaBarFromXiBarPx=
 860   fLambdaBarFromXiBar->ProjectionX("Bar"); lambdaBarFromXiBarPx->Sumw2();
 861
 862   fLambdaBarMC = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaBarMC")) ;
 863   if (!fLambdaBarMC) {
 864      Printf("ERROR: fLambdaBarMC not available");
 865      return;
 866   }
 867   TH1D *lambdaBarMcPx=fLambdaBarMC->ProjectionX(); lambdaBarMcPx->Sumw2();
 868
 869   fLambdaBarAs = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaBarAs")) ;
 870   if (!fLambdaBarAs) {
 871      Printf("ERROR: fLambdaBarAs not available");
 872      return;
 873   }
 874   TH1D *lambdaBarAsPx=fLambdaBarAs->ProjectionX();
 875   lambdaBarAsPx->Sumw2(); //lambdaBarAsPx->Scale(0.64);
 876
 877   fLambdaBarSi = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaBarSi")) ;
 878   if (!fLambdaBarSi) {
 879      Printf("ERROR: fLambdaBarSi not available");
 880      return;
 881   }
 882   TH1D *lambdaBarSiPx=fLambdaBarSi->ProjectionX();
 883   lambdaBarSiPx->SetName("fLambdaBarPt");
 884   lambdaBarSiPx->Sumw2();
 885
 886   fLambdaBarEff = dynamic_cast<TH1D*>(fOutput->FindObject("fLambdaBarEff")) ;
 887   if (!fLambdaBarEff) {
 888      Printf("ERROR: fLambdaBarEff not available");
 889      return;
 890   }
 891   fLambdaBarPt = dynamic_cast<TH1D*>(fOutput->FindObject("fLambdaBarPt")) ;
 892   if (!fLambdaBarPt) {
 893      Printf("ERROR: fLambdaBarPt not available");
 894      return;
 895   }
 896
 897
 898   if (!gROOT->IsBatch()) {
 899
 900     TCanvas *c1 = new TCanvas("c1","Mulitplicity");
 901     c1->SetLogy();
 902     fMult->DrawCopy() ;
 903
 904     new TCanvas("c2","dE/dx");
 905     fdEdx->DrawCopy() ;
 906
 907     new TCanvas("c3","dE/dx with PID");
 908     fdEdxPid->DrawCopy() ;
 909
 910     if (fIsMC) {
 911        /*
 912        TH1D effK(*k0sAsPx); effK.SetTitle("Efficiency for K0s");
 913        effK.Divide(k0sAsPx,k0sMcPx,1,1,"b");
 914        new TCanvas("c4","Efficiency for K0s");
 915        effK.DrawCopy("E") ;
 916        */
 917
 918        //+++ Lambdas
 919        fLambdaEff->Divide(lambdaAsPx,lambdaMcPx,1,1,"b");
 920        new TCanvas("c5","Efficiency for #Lambda");
 921        fLambdaEff->DrawCopy("E") ;
 922
 923        lambdaSiPx->Add(lambdaFromXiPx,-1);
 924        lambdaSiPx->Divide(fLambdaEff);
 925
 926        new TCanvas("c6","Corrected #Lambda pt");
 927        lambdaSiPx->SetTitle("Corrected #Lambda pt");
 928       *fLambdaPt = *lambdaSiPx;
 929        fLambdaPt->SetLineColor(2);
 930        fLambdaPt->DrawCopy("E");
 931
 932        lambdaMcPx->DrawCopy("same");
 933
 934
 935        //+++ anti-Lambdas
 936        fLambdaBarEff->Divide(lambdaBarAsPx,lambdaBarMcPx,1,1,"b");
 937        new TCanvas("c7","Efficiency for anti-#Lambda");
 938        fLambdaBarEff->DrawCopy("E") ;
 939
 940        lambdaBarSiPx->Add(lambdaBarFromXiBarPx,-1);
 941        lambdaBarSiPx->Divide(fLambdaBarEff);
 942
 943        new TCanvas("c8","Corrected anti-#Lambda pt");
 944        lambdaBarSiPx->SetTitle("Corrected anti-#Lambda pt");
 945       *fLambdaBarPt = *lambdaBarSiPx;
 946        fLambdaBarPt->SetLineColor(2);
 947        fLambdaBarPt->DrawCopy("E");
 948
 949        lambdaBarMcPx->DrawCopy("same");
 950     } else {
 951        new TCanvas("c6","Raw #Lambda pt");
 952        lambdaSiPx->SetTitle("Raw #Lambda pt");
 953       *fLambdaPt = *lambdaSiPx;
 954        fLambdaPt->SetLineColor(2);
 955        fLambdaPt->DrawCopy("E");
 956
 957
 958        new TCanvas("c7","Raw anti-#Lambda pt");
 959        lambdaBarSiPx->SetTitle("Raw anti-#Lambda pt");
 960       *fLambdaBarPt = *lambdaBarSiPx;
 961        fLambdaBarPt->SetLineColor(2);
 962        fLambdaBarPt->DrawCopy("E");
 963     }
 964   }
 965 }