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