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","L_{T} vs p_{T} 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","L_{T} vs p_{T} 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","L_{T} vs p_{T} 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","L_{T} vs p_{T} 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
 276   Double_t cpa=v0->GetV0CosineOfPointingAngle();
 277   if (cpa<0.998) return kFALSE;
 278
 279   Double_t xx,yy,zz; v0->GetXYZ(xx,yy,zz);
 280   Double_t r2=xx*xx + yy*yy;
 281   if (r2<0.9*0.9) return kFALSE;
 282   if (r2>100*100) return kFALSE;
 283
 284   return kTRUE;
 285 }
 286
 287 static Bool_t AcceptCascade(const AliESDcascade *cs, const AliESDEvent *esd) {
 288
 289   Int_t bidx=TMath::Abs(cs->GetBindex());
 290   AliESDtrack *btrack=esd->GetTrack(bidx);
 291   if (!AcceptTrack(btrack)) return kFALSE;
 292
 293   Float_t xy,z0;
 294   btrack->GetImpactParameters(xy,z0);
 295   if (TMath::Abs(xy)<0.03) return kFALSE;
 296
 297   const AliESDVertex *vtx=esd->GetPrimaryVertexSPD();
 298   if (!vtx->GetStatus()) {
 299      vtx=esd->GetPrimaryVertexTracks();
 300      if (!vtx->GetStatus()) return kFALSE;
 301   }
 302   Double_t xv=vtx->GetXv(), yv=vtx->GetYv(), zv=vtx->GetZv();
 303   if (cs->GetCascadeCosineOfPointingAngle(xv,yv,zv) < 0.999) return kFALSE;
 304
 305   if (cs->GetDcaXiDaughters() > 0.3) return kFALSE;
 306
 307   return kTRUE;
 308 }
 309
 310 static Bool_t AcceptPID(const AliPIDResponse *pidResponse,
 311                         const AliESDtrack *ptrack, AliStack *stack) {
 312
 313   const AliExternalTrackParam *par=ptrack->GetInnerParam();
 314   if (!par) return kTRUE;
 315   if (par->GetP() > 1.) return kTRUE;
 316
 317
 318   if (stack) {
 319     // MC PID
 320     Int_t plab=TMath::Abs(ptrack->GetLabel());
 321     TParticle *pp=stack->Particle(plab);
 322     if (!pp) return kTRUE;
 323     if (pp->GetPDG()->Charge() > 0) {
 324        if (pp->GetPdgCode() == kProton)    return kTRUE;
 325     } else {
 326        if (pp->GetPdgCode() == kProtonBar) return kTRUE;
 327     }
 328   } else {
 329     // Real PID
 330     Double_t nsig=pidResponse->NumberOfSigmasTPC(ptrack,AliPID::kProton);
 331     if (TMath::Abs(nsig) < 3.) return kTRUE;
 332   }
 333
 334   return kFALSE;
 335 }
 336
 337 static TParticle*
 338 AssociateV0(const AliESDtrack *ptrack,const AliESDtrack *ntrack,AliStack *stack,
 339 TParticle *&mcp) {
 340   //
 341   // Try to associate the V0 with the daughters ptrack and ntrack
 342   // with the Monte Carlo
 343   //
 344   if (!stack) return 0;
 345
 346   Int_t nlab=TMath::Abs(ntrack->GetLabel());
 347   TParticle *n=stack->Particle(nlab);
 348   if (!n) return 0;
 349
 350   Int_t plab=TMath::Abs(ptrack->GetLabel());
 351   TParticle *p=stack->Particle(plab);
 352   if (!p) return 0;
 353
 354   Int_t imp=p->GetFirstMother();
 355   if (imp<0) return 0;
 356   if (imp>=stack->GetNtrack()) return 0;
 357   TParticle *p0=stack->Particle(imp); // V0 particle, mother of the pos. track
 358   if (!p0) return 0;
 359
 360   Int_t imn=n->GetFirstMother();
 361   if (imp != imn) {  // Check decays of the daughters
 362      return 0; // Fixme
 363   }
 364
 365   Int_t code=p0->GetPdgCode();
 366   if (code != kK0Short)
 367      if (code != kLambda0)
 368         if (code != kLambda0Bar) return 0;
 369
 370   mcp=p;
 371
 372   return p0;
 373 }
 374
 375
 376 void AliAnalysisTaskCTauPbPb::UserExec(Option_t *)
 377 {
 378   const Double_t yMax=0.5;
 379   const Double_t pMin=0.0;
 380   const Double_t lMax=0.001;
 381
 382   AliESDEvent *esd=(AliESDEvent *)InputEvent();
 383
 384   if (!esd) {
 385     Printf("ERROR: esd not available");
 386     return;
 387   }
 388
 389   // Vertex selection
 390   const AliESDVertex *vtx=esd->GetPrimaryVertex();
 391   if (!vtx->GetStatus()) return;
 392   Double_t xv=vtx->GetXv(), yv=vtx->GetYv(), zv=vtx->GetZv();
 393
 394   if (TMath::Abs(zv) > 10.) return ;
 395
 396
 397   // Physics selection
 398   AliAnalysisManager *mgr= AliAnalysisManager::GetAnalysisManager();
 399   AliInputEventHandler *hdr=(AliInputEventHandler*)mgr->GetInputEventHandler();
 400   UInt_t maskIsSelected = hdr->IsEventSelected();
 401   Bool_t isSelected = (maskIsSelected & AliVEvent::kMB);
 402   if (!isSelected) return;
 403
 404
 405   fMult->Fill(-100); //event counter
 406
 407
 408   // Centrality selection
 409   AliCentrality *cent=esd->GetCentrality();
 410   if (!cent->IsEventInCentralityClass(fCMin,fCMax,"V0M")) return;
 411
 412   AliPIDResponse *pidResponse = hdr->GetPIDResponse();
 413
 414
 415   //+++++++ MC
 416   AliStack *stack = 0x0;
 417   Double_t mcXv=0., mcYv=0., mcZv=0.;
 418
 419   if (fIsMC) {
 420      AliMCEvent *mcEvent = MCEvent();
 421      stack = mcEvent->Stack();
 422      if (!stack) {
 423         Printf("ERROR: stack not available");
 424         return;
 425      }
 426
 427      const AliVVertex *mcVtx=mcEvent->GetPrimaryVertex();
 428
 429      mcXv=mcVtx->GetX(); mcYv=mcVtx->GetY(); mcZv=mcVtx->GetZ();
 430
 431      Int_t ntrk=stack->GetNtrack();
 432      while (ntrk--) {
 433        TParticle *p0=stack->Particle(ntrk);
 434        Int_t code=p0->GetPdgCode();
 435        if (code != kK0Short)
 436          if (code != kLambda0)
 437             if (code != kLambda0Bar) continue;
 438
 439        Int_t plab=p0->GetFirstDaughter(), nlab=p0->GetLastDaughter();
 440        if (nlab==plab) continue;
 441        TParticle *part = stack->Particle(plab);
 442        if (!part) continue;
 443        TParticlePDG *partPDG = part->GetPDG();
 444        if (!partPDG) continue;
 445        Double_t charge=partPDG->Charge();
 446        if (charge == 0.) continue;
 447
 448        Double_t pt=p0->Pt();
 449        if (pt<pMin) continue;
 450        if (TMath::Abs(p0->Y())>yMax) continue;
 451
 452        Double_t x=p0->Vx(), y=p0->Vy(), z=p0->Vz();
 453        Double_t dx=mcXv-x, dy=mcYv-y, dz=mcZv-z;
 454        Double_t l=TMath::Sqrt(dx*dx + dy*dy + dz*dz);
 455
 456        if (l > lMax) continue; // secondary V0
 457
 458        x=part->Vx(); y=part->Vy();
 459        dx=mcXv-x; dy=mcYv-y;
 460        Double_t lt=TMath::Sqrt(dx*dx + dy*dy);
 461
 462        switch (code) {
 463        case kK0Short:
 464           fK0sMC->Fill(pt,lt);
 465           break;
 466        case kLambda0:
 467           fLambdaMC->Fill(pt,lt);
 468           break;
 469        case kLambda0Bar:
 470           fLambdaBarMC->Fill(pt,lt);
 471           break;
 472        default: break;
 473        }
 474      }
 475   }
 476   //+++++++
 477
 478
 479   Int_t ntrk1=esd->GetNumberOfTracks();
 480   Int_t mult=0;
 481   for (Int_t i=0; i<ntrk1; i++) {
 482     AliESDtrack *t=esd->GetTrack(i);
 483     if (!t->IsOn(AliESDtrack::kTPCrefit)) continue;
 484     Float_t xy,z0;
 485     t->GetImpactParameters(xy,z0);
 486     if (TMath::Abs(xy)>3.) continue;
 487     if (TMath::Abs(z0)>3.) continue;
 488     Double_t pt=t->Pt(),pz=t->Pz();
 489     if (TMath::Abs(pz/pt)>0.8) continue;
 490     mult++;
 491
 492     Double_t p=t->GetInnerParam()->GetP();
 493     Double_t dedx=t->GetTPCsignal()/47.;
 494     fdEdx->Fill(p,dedx,1);
 495
 496     Double_t nsig=pidResponse->NumberOfSigmasTPC(t,AliPID::kProton);
 497     if (TMath::Abs(nsig) < 3.) fdEdxPid->Fill(p,dedx,1);
 498
 499   }
 500   fMult->Fill(mult);
 501
 502
 503   Int_t nv0 = esd->GetNumberOfV0s();
 504   while (nv0--) {
 505       AliESDv0 *v0=esd->GetV0(nv0);
 506
 507       Double_t pt=v0->Pt();
 508       if (pt < pMin) continue;
 509
 510       if (!AcceptV0(v0,esd)) continue;
 511
 512       Int_t nidx=TMath::Abs(v0->GetNindex());
 513       AliESDtrack *ntrack=esd->GetTrack(nidx);
 514       Int_t pidx=TMath::Abs(v0->GetPindex());
 515       AliESDtrack *ptrack=esd->GetTrack(pidx);
 516
 517       Double_t x,y,z; v0->GetXYZ(x,y,z);
 518       Double_t dx1=x-xv, dy1=y-yv;
 519       Double_t lt=TMath::Sqrt(dx1*dx1 + dy1*dy1);
 520
 521       if (lt/pt > 3*7.89/1.1157) continue;
 522
 523       //--- V0 switches
 524       Bool_t isK0s=kTRUE;
 525       Bool_t isLambda=kTRUE;
 526       Bool_t isLambdaBar=kTRUE;
 527
 528       if (0.4977*lt/pt > 3*2.68) isK0s=kFALSE;
 529       if (1.1157*lt/pt > 3*7.89) isLambdaBar=isLambda=kFALSE;
 530
 531       if (v0->PtArmV0() < 0.2*TMath::Abs(v0->AlphaV0())) isK0s=kFALSE;
 532
 533       if (!AcceptPID(pidResponse, ptrack, stack)) isLambda=kFALSE;
 534       if (!AcceptPID(pidResponse, ntrack, stack)) isLambdaBar=kFALSE;
 535
 536       Double_t yK0s=TMath::Abs(v0->RapK0Short());
 537       Double_t yLam=TMath::Abs(v0->RapLambda());
 538       if (yK0s > yMax) isK0s=kFALSE;
 539       if (yLam > yMax) isLambda=isLambdaBar=kFALSE;
 540       //---
 541
 542       Double_t mass=0., m=0., s=0.;
 543       if (isK0s) {
 544          v0->ChangeMassHypothesis(kK0Short);
 545
 546          mass=v0->GetEffMass();
 547          fK0sM->Fill(mass,pt);
 548
 549          m=TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass();
 550          s=0.0044 + (0.008-0.0044)/(10-1)*(pt - 1.);
 551          if (TMath::Abs(m-mass) < 3*s) {
 552             fK0sSi->Fill(pt,lt);
 553          } else {
 554             isK0s=kFALSE;
 555          }
 556          if (TMath::Abs(m-mass + 4.5*s) < 1.5*s) {
 557             fK0sSi->Fill(pt,lt,-1);
 558          }
 559          if (TMath::Abs(m-mass - 4.5*s) < 1.5*s) {
 560             fK0sSi->Fill(pt,lt,-1);
 561          }
 562       }
 563
 564       if (isLambda) {
 565          v0->ChangeMassHypothesis(kLambda0);
 566
 567          mass=v0->GetEffMass();
 568          fLambdaM->Fill(mass,pt);
 569
 570          m=TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass();
 571          s=0.0023 + (0.004-0.0023)/(6-1)*(pt-1);
 572          if (TMath::Abs(m-mass) < 3*s) {
 573             fLambdaSi->Fill(pt,lt);
 574          } else {
 575             isLambda=kFALSE;
 576          }
 577          if (TMath::Abs(m-mass + 4.5*s) < 1.5*s) {
 578             fLambdaSi->Fill(pt,lt,-1);
 579          }
 580          if (TMath::Abs(m-mass - 4.5*s) < 1.5*s) {
 581             fLambdaSi->Fill(pt,lt,-1);
 582          }
 583       }
 584
 585       if (isLambdaBar) {
 586          v0->ChangeMassHypothesis(kLambda0Bar);
 587
 588          mass=v0->GetEffMass();
 589          fLambdaBarM->Fill(mass,pt);
 590
 591          m=TDatabasePDG::Instance()->GetParticle(kLambda0Bar)->Mass();
 592          s=0.0023 + (0.004-0.0023)/(6-1)*(pt-1);
 593          if (TMath::Abs(m-mass) < 3*s) {
 594             fLambdaBarSi->Fill(pt,lt);
 595          } else {
 596             isLambdaBar=kFALSE;
 597          }
 598          if (TMath::Abs(m-mass + 4.5*s) < 1.5*s) {
 599             fLambdaBarSi->Fill(pt,lt,-1);
 600          }
 601          if (TMath::Abs(m-mass - 4.5*s) < 1.5*s) {
 602             fLambdaBarSi->Fill(pt,lt,-1);
 603          }
 604       }
 605
 606       if (!fIsMC) continue;
 607
 608       //++++++ MC
 609       if (!isK0s)
 610          if (!isLambda)
 611              if (!isLambdaBar) continue;//check MC only for the accepted V0s
 612
 613       TParticle *mcp=0;
 614       TParticle *mc0=AssociateV0(ptrack,ntrack,stack,mcp);
 615       if (!mc0) continue;
 616
 617       Double_t ptAs=mc0->Pt();
 618       if (ptAs < pMin) continue;
 619       Double_t yAs=mc0->Y();
 620       if (TMath::Abs(yAs) > yMax) continue;
 621
 622       Int_t code=mc0->GetPdgCode();
 623
 624       Double_t dx = mcXv - mcp->Vx(), dy = mcYv - mcp->Vy();
 625       Double_t ltAs=TMath::Sqrt(dx*dx + dy*dy);
 626
 627       Double_t dz=mcZv - mc0->Vz(); dy=mcYv - mc0->Vy(); dx=mcXv - mc0->Vx();
 628       Double_t l = TMath::Sqrt(dx*dx + dy*dy + dz*dz);
 629       if (l<lMax) { // Primary V0s
 630          switch (code) {
 631          case kK0Short:
 632             if (isK0s)       fK0sAs->Fill(ptAs,ltAs);
 633             break;
 634          case kLambda0:
 635             if (isLambda)    fLambdaAs->Fill(ptAs,ltAs);
 636             break;
 637          case kLambda0Bar:
 638             if (isLambdaBar) fLambdaBarAs->Fill(ptAs,ltAs);
 639             break;
 640          default: break;
 641          }
 642       } else {
 643          if (code==kK0Short) continue;
 644
 645          Int_t nx=mc0->GetFirstMother();
 646          if (nx<0) continue;
 647          if (nx>=stack->GetNtrack()) continue;
 648          TParticle *xi=stack->Particle(nx);
 649          if (!xi) continue;
 650          Int_t xcode=xi->GetPdgCode();
 651
 652          switch (code) {
 653          case kLambda0:
 654             if (isLambda)
 655             if ((xcode==kXiMinus) || (xcode==3322))
 656                fLambdaFromXi->Fill(ptAs,ltAs,xi->Pt());
 657             break;
 658          case kLambda0Bar:
 659             if (isLambdaBar)
 660             if ((xcode==kXiPlusBar)||(xcode==-3322))
 661                fLambdaBarFromXiBar->Fill(ptAs,ltAs,xi->Pt());
 662             break;
 663          default: break;
 664          }
 665       }
 666       //++++++
 667
 668   }
 669
 670   Double_t kine0;
 671   Int_t ncs=esd->GetNumberOfCascades();
 672   for (Int_t i=0; i<ncs; i++) {
 673       AliESDcascade *cs=esd->GetCascade(i);
 674
 675       Double_t pt=cs->Pt();
 676       if (pt < pMin) continue;
 677       if (TMath::Abs(cs->RapXi()) > yMax) continue;
 678       if (!AcceptCascade(cs,esd)) continue;
 679
 680       AliESDv0 *v0 = (AliESDv0*)cs;
 681       if (v0->Pt() < pMin) continue;
 682       if (TMath::Abs(v0->RapLambda()) > yMax) continue;
 683       if (!AcceptV0(v0,esd)) continue;
 684
 685       //--- Cascade switches
 686       Bool_t isXiMinus=kTRUE;
 687       Bool_t isXiPlusBar=kTRUE;
 688
 689       Int_t pidx=TMath::Abs(v0->GetPindex());
 690       AliESDtrack *ptrack=esd->GetTrack(pidx);
 691       if (!AcceptPID(pidResponse, ptrack, stack)) isXiMinus=kFALSE;
 692
 693       Int_t nidx=TMath::Abs(v0->GetNindex());
 694       AliESDtrack *ntrack=esd->GetTrack(nidx);
 695       if (!AcceptPID(pidResponse, ntrack, stack)) isXiPlusBar=kFALSE;
 696
 697       Int_t charge=cs->Charge();
 698       if (charge > 0) isXiMinus=kFALSE;
 699       if (charge < 0) isXiPlusBar=kFALSE;
 700       //---
 701
 702       if (isXiMinus) {
 703          cs->ChangeMassHypothesis(kine0,kXiMinus);
 704          Double_t mass=cs->GetEffMassXi();
 705          fXiM->Fill(mass,pt);
 706          Double_t m=TDatabasePDG::Instance()->GetParticle(kXiMinus)->Mass();
 707          //Double_t s=0.0037;
 708          Double_t s=0.002 + (0.0032-0.002)/(6-1.5)*(pt-1.5);
 709          if (TMath::Abs(m-mass) < 3*s) {
 710             fXiSiP->Fill(pt);
 711          } else {
 712             isXiMinus=kFALSE;
 713          }
 714          if (TMath::Abs(m-mass + 4.5*s) < 1.5*s) {
 715             fXiSiP->Fill(pt,-1);
 716          }
 717          if (TMath::Abs(m-mass - 4.5*s) < 1.5*s) {
 718             fXiSiP->Fill(pt,-1);
 719          }
 720       }
 721
 722       if (isXiPlusBar) {
 723          cs->ChangeMassHypothesis(kine0,kXiPlusBar);
 724          Double_t mass=cs->GetEffMassXi();
 725          fXiBarM->Fill(mass,pt);
 726          Double_t m=TDatabasePDG::Instance()->GetParticle(kXiPlusBar)->Mass();
 727          //Double_t s=0.0037;
 728          Double_t s=0.002 + (0.0032-0.002)/(6-1.5)*(pt-1.5);
 729          if (TMath::Abs(m-mass) < 3*s) {
 730             fXiBarSiP->Fill(pt);
 731          } else {
 732             isXiPlusBar=kFALSE;
 733          }
 734          if (TMath::Abs(m-mass + 4.5*s) < 1.5*s) {
 735             fXiBarSiP->Fill(pt,-1);
 736          }
 737          if (TMath::Abs(m-mass - 4.5*s) < 1.5*s) {
 738             fXiBarSiP->Fill(pt,-1);
 739          }
 740       }
 741
 742       if (!fIsMC) continue;
 743
 744       //++++++ MC
 745       if (!isXiMinus)
 746          if (!isXiPlusBar) continue;//check MC only for the accepted cascades
 747       // Here is the future association with MC
 748   }
 749
 750 }
 751
 752 void AliAnalysisTaskCTauPbPb::Terminate(Option_t *)
 753 {
 754    // The Terminate() function is the last function to be called during
 755    // a query. It always runs on the client, it can be used to present
 756    // the results graphically or save the results to file.
 757
 758   fOutput=(TList*)GetOutputData(1);
 759   if (!fOutput) {
 760      Printf("ERROR: fOutput not available");
 761      return;
 762   }
 763
 764   fMult = dynamic_cast<TH1F*>(fOutput->FindObject("fMult")) ;
 765   if (!fMult) {
 766      Printf("ERROR: fMult not available");
 767      return;
 768   }
 769
 770   fdEdx = dynamic_cast<TH2F*>(fOutput->FindObject("fdEdx")) ;
 771   if (!fdEdx) {
 772      Printf("ERROR: fdEdx not available");
 773      return;
 774   }
 775
 776   fdEdxPid = dynamic_cast<TH2F*>(fOutput->FindObject("fdEdxPid")) ;
 777   if (!fdEdxPid) {
 778      Printf("ERROR: fdEdxPid not available");
 779      return;
 780   }
 781
 782
 783   fK0sMC = dynamic_cast<TH2F*>(fOutput->FindObject("fK0sMC")) ;
 784   if (!fK0sMC) {
 785      Printf("ERROR: fK0sMC not available");
 786      return;
 787   }
 788   TH1D *k0sMcPx=fK0sMC->ProjectionX(); k0sMcPx->Sumw2();
 789   fK0sAs = dynamic_cast<TH2F*>(fOutput->FindObject("fK0sAs")) ;
 790   if (!fK0sAs) {
 791      Printf("ERROR: fK0sAs not available");
 792      return;
 793   }
 794   TH1D *k0sAsPx=fK0sAs->ProjectionX();
 795   k0sAsPx->Sumw2(); //k0sAsPx->Scale(0.69);
 796
 797
 798
 799   // Lambda histograms
 800   fLambdaFromXi = dynamic_cast<TH3F*>(fOutput->FindObject("fLambdaFromXi")) ;
 801   if (!fLambdaFromXi) {
 802      Printf("ERROR: fLambdaFromXi not available");
 803      return;
 804   }
 805   TH1D *lambdaFromXiPx=fLambdaFromXi->ProjectionX(); lambdaFromXiPx->Sumw2();
 806
 807   fLambdaMC = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaMC")) ;
 808   if (!fLambdaMC) {
 809      Printf("ERROR: fLambdaMC not available");
 810      return;
 811   }
 812   TH1D *lambdaMcPx=fLambdaMC->ProjectionX(); lambdaMcPx->Sumw2();
 813
 814   fLambdaAs = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaAs")) ;
 815   if (!fLambdaAs) {
 816      Printf("ERROR: fLambdaAs not available");
 817      return;
 818   }
 819   TH1D *lambdaAsPx=fLambdaAs->ProjectionX();
 820   lambdaAsPx->Sumw2(); //lambdaAsPx->Scale(0.64);
 821
 822   fLambdaSi = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaSi")) ;
 823   if (!fLambdaSi) {
 824      Printf("ERROR: fLambdaSi not available");
 825      return;
 826   }
 827   TH1D *lambdaSiPx=fLambdaSi->ProjectionX();
 828   lambdaSiPx->SetName("fLambdaPt");
 829   lambdaSiPx->Sumw2();
 830
 831   fLambdaEff = dynamic_cast<TH1D*>(fOutput->FindObject("fLambdaEff")) ;
 832   if (!fLambdaEff) {
 833      Printf("ERROR: fLambdaEff not available");
 834      return;
 835   }
 836   fLambdaPt = dynamic_cast<TH1D*>(fOutput->FindObject("fLambdaPt")) ;
 837   if (!fLambdaPt) {
 838      Printf("ERROR: fLambdaPt not available");
 839      return;
 840   }
 841
 842
 843   // anti-Lambda histograms
 844   fLambdaBarFromXiBar =
 845   dynamic_cast<TH3F*>(fOutput->FindObject("fLambdaBarFromXiBar")) ;
 846   if (!fLambdaBarFromXiBar) {
 847      Printf("ERROR: fLambdaBarFromXiBar not available");
 848      return;
 849   }
 850   TH1D *lambdaBarFromXiBarPx=
 851   fLambdaBarFromXiBar->ProjectionX("Bar"); lambdaBarFromXiBarPx->Sumw2();
 852
 853   fLambdaBarMC = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaBarMC")) ;
 854   if (!fLambdaBarMC) {
 855      Printf("ERROR: fLambdaBarMC not available");
 856      return;
 857   }
 858   TH1D *lambdaBarMcPx=fLambdaBarMC->ProjectionX(); lambdaBarMcPx->Sumw2();
 859
 860   fLambdaBarAs = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaBarAs")) ;
 861   if (!fLambdaBarAs) {
 862      Printf("ERROR: fLambdaBarAs not available");
 863      return;
 864   }
 865   TH1D *lambdaBarAsPx=fLambdaBarAs->ProjectionX();
 866   lambdaBarAsPx->Sumw2(); //lambdaBarAsPx->Scale(0.64);
 867
 868   fLambdaBarSi = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaBarSi")) ;
 869   if (!fLambdaBarSi) {
 870      Printf("ERROR: fLambdaBarSi not available");
 871      return;
 872   }
 873   TH1D *lambdaBarSiPx=fLambdaBarSi->ProjectionX();
 874   lambdaBarSiPx->SetName("fLambdaBarPt");
 875   lambdaBarSiPx->Sumw2();
 876
 877   fLambdaBarEff = dynamic_cast<TH1D*>(fOutput->FindObject("fLambdaBarEff")) ;
 878   if (!fLambdaBarEff) {
 879      Printf("ERROR: fLambdaBarEff not available");
 880      return;
 881   }
 882   fLambdaBarPt = dynamic_cast<TH1D*>(fOutput->FindObject("fLambdaBarPt")) ;
 883   if (!fLambdaBarPt) {
 884      Printf("ERROR: fLambdaBarPt not available");
 885      return;
 886   }
 887
 888
 889   if (!gROOT->IsBatch()) {
 890
 891     TCanvas *c1 = new TCanvas("c1","Mulitplicity");
 892     c1->SetLogy();
 893     fMult->DrawCopy() ;
 894
 895     new TCanvas("c2","dE/dx");
 896     fdEdx->DrawCopy() ;
 897
 898     new TCanvas("c3","dE/dx with PID");
 899     fdEdxPid->DrawCopy() ;
 900
 901     if (fIsMC) {
 902        /*
 903        TH1D effK(*k0sAsPx); effK.SetTitle("Efficiency for K0s");
 904        effK.Divide(k0sAsPx,k0sMcPx,1,1,"b");
 905        new TCanvas("c4","Efficiency for K0s");
 906        effK.DrawCopy("E") ;
 907        */
 908
 909        //+++ Lambdas
 910        fLambdaEff->Divide(lambdaAsPx,lambdaMcPx,1,1,"b");
 911        new TCanvas("c5","Efficiency for #Lambda");
 912        fLambdaEff->DrawCopy("E") ;
 913
 914        lambdaSiPx->Add(lambdaFromXiPx,-1);
 915        lambdaSiPx->Divide(fLambdaEff);
 916
 917        new TCanvas("c6","Corrected #Lambda pt");
 918        lambdaSiPx->SetTitle("Corrected #Lambda pt");
 919       *fLambdaPt = *lambdaSiPx;
 920        fLambdaPt->SetLineColor(2);
 921        fLambdaPt->DrawCopy("E");
 922
 923        lambdaMcPx->DrawCopy("same");
 924
 925
 926        //+++ anti-Lambdas
 927        fLambdaBarEff->Divide(lambdaBarAsPx,lambdaBarMcPx,1,1,"b");
 928        new TCanvas("c7","Efficiency for anti-#Lambda");
 929        fLambdaBarEff->DrawCopy("E") ;
 930
 931        lambdaBarSiPx->Add(lambdaBarFromXiBarPx,-1);
 932        lambdaBarSiPx->Divide(fLambdaBarEff);
 933
 934        new TCanvas("c8","Corrected anti-#Lambda pt");
 935        lambdaBarSiPx->SetTitle("Corrected anti-#Lambda pt");
 936       *fLambdaBarPt = *lambdaBarSiPx;
 937        fLambdaBarPt->SetLineColor(2);
 938        fLambdaBarPt->DrawCopy("E");
 939
 940        lambdaBarMcPx->DrawCopy("same");
 941     } else {
 942        new TCanvas("c6","Raw #Lambda pt");
 943        lambdaSiPx->SetTitle("Raw #Lambda pt");
 944       *fLambdaPt = *lambdaSiPx;
 945        fLambdaPt->SetLineColor(2);
 946        fLambdaPt->DrawCopy("E");
 947
 948
 949        new TCanvas("c7","Raw anti-#Lambda pt");
 950        lambdaBarSiPx->SetTitle("Raw anti-#Lambda pt");
 951       *fLambdaBarPt = *lambdaBarSiPx;
 952        fLambdaBarPt->SetLineColor(2);
 953        fLambdaBarPt->DrawCopy("E");
 954     }
 955   }
 956 }