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 fOutput(0),
  47 fMult(0),
  48 fCosPA(0),
  49 fDtrDCA(0),
  50 fdEdx(0),
  51 fdEdxPid(0),
  52
  53 fK0sM(0),
  54 fK0sSi(0),
  55 fK0sMC(0),
  56 fK0sAs(0),
  57
  58 fLambdaM(0),
  59 fLambdaSi(0),
  60 fLambdaMC(0),
  61 fLambdaAs(0),
  62
  63 fLambdaEff(0),
  64 fLambdaPt(0),
  65
  66 fLambdaBarM(0),
  67 fLambdaBarSi(0),
  68 fLambdaBarMC(0),
  69 fLambdaBarAs(0),
  70
  71 fLambdaBarEff(0),
  72 fLambdaBarPt(0),
  73
  74 fLambdaFromXi(0),
  75 fXiM(0),
  76 fXiSiP(0),
  77
  78 fLambdaBarFromXiBar(0),
  79 fXiBarM(0),
  80 fXiBarSiP(0)
  81 {
  82   // Constructor. Initialization of pointers
  83   DefineOutput(1, TList::Class());
  84 }
  85
  86 void AliAnalysisTaskCTauPbPbaod::UserCreateOutputObjects()
  87 {
  88   Int_t    nbins=100;  // number of bins
  89   Double_t ltMax=100.;
  90   Double_t ptMax=10.;
  91
  92   fOutput = new TList();
  93   fOutput->SetOwner();
  94
  95
  96   fMult=new TH1F("fMult","Multiplicity",1100,0.,3300);
  97   fMult->GetXaxis()->SetTitle("N tracks");
  98   fOutput->Add(fMult);
  99
 100   fCosPA=new TH1F("fCosPA","Cos(PA) distribution",50,0.9975,1.0005);
 101   fCosPA->GetXaxis()->SetTitle("Cos(PA)");
 102   fOutput->Add(fCosPA);
 103
 104   fDtrDCA=new TH1F("fDtrDCA","DCA between V0 daughters",50,0.0,1.5);
 105   fDtrDCA->GetXaxis()->SetTitle("DCA (rel. u.)");
 106   fOutput->Add(fDtrDCA);
 107
 108   fdEdx=new TH2F("fdEdx","dE/dx",50,0.2,3,50,0.,6.);
 109   fOutput->Add(fdEdx);
 110
 111   fdEdxPid=new TH2F("fdEdxPid","dE/dx with PID",50,0.2,3,50,0.,6.);
 112   fOutput->Add(fdEdxPid);
 113
 114   fK0sM =
 115   new TH2F("fK0sM", "Mass for K^{0}_{s}", nbins/2,0.448,0.548,nbins,0.,ptMax);
 116   fK0sM->GetXaxis()->SetTitle("Mass (GeV/c)");
 117   fOutput->Add(fK0sM);
 118
 119   fK0sSi =
 120   new TH2F("fK0sSi","L_{T} vs p_{T} for K^{0}_{s}, side-band subtracted",
 121   nbins,0.,ptMax,nbins,0.,ltMax);
 122   fK0sSi->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 123   fK0sSi->GetYaxis()->SetTitle("L_{T} (cm)");
 124   fOutput->Add(fK0sSi);
 125
 126   fK0sMC =
 127   new TH2F("fK0sMC","L_{T} vs p_{T} for K^{0}_{s}, from MC stack",
 128   nbins,0.,ptMax,nbins,0.,ltMax);
 129   fK0sMC->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 130   fK0sMC->GetYaxis()->SetTitle("L_{T} (cm)");
 131   fOutput->Add(fK0sMC);
 132
 133   fK0sAs =
 134   new TH2F("fK0sAs", "L_{T} vs p_{T} for K^{0}_{s}, associated",
 135   nbins,0.,ptMax,nbins,0.,ltMax);
 136   fK0sAs->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 137   fK0sAs->GetYaxis()->SetTitle("L_{T} (cm)");
 138   fOutput->Add(fK0sAs);
 139
 140   //----------------------
 141
 142   fLambdaM =
 143   new TH2F("fLambdaM","Mass for \\Lambda", nbins, 1.065, 1.165,nbins,0.,ptMax);
 144   fLambdaM->GetXaxis()->SetTitle("Mass (GeV/c)");
 145   fOutput->Add(fLambdaM);
 146
 147   fLambdaSi =
 148   new TH2F("fLambdaSi","L_{T} vs p_{T} for \\Lambda, side-band subtructed",
 149   nbins,0.,ptMax,nbins,0.,ltMax);
 150   fLambdaSi->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 151   fLambdaSi->GetYaxis()->SetTitle("L_{T} (cm)");
 152   fOutput->Add(fLambdaSi);
 153
 154   fLambdaMC =
 155   new TH2F("fLambdaMC","L_{T} vs p_{T} for \\Lambda, from MC stack",
 156   nbins,0.,ptMax,nbins,0.,ltMax);
 157   fLambdaMC->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 158   fLambdaMC->GetYaxis()->SetTitle("L_{T} (cm)");
 159   fOutput->Add(fLambdaMC);
 160
 161   fLambdaAs =
 162   new TH2F("fLambdaAs","L_{T} vs p_{T} for \\Lambda, associated",
 163   nbins,0.,ptMax,nbins,0.,ltMax);
 164   fLambdaAs->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 165   fLambdaAs->GetYaxis()->SetTitle("L_{T} (cm)");
 166   fOutput->Add(fLambdaAs);
 167
 168   //----------------------
 169
 170   fLambdaEff=fLambdaAs->ProjectionX();
 171   fLambdaEff->SetName("fLambdaEff");
 172   fLambdaEff->SetTitle("Efficiency for #Lambda");
 173   fOutput->Add(fLambdaEff);
 174
 175   fLambdaPt=fLambdaAs->ProjectionX();
 176   fLambdaPt->SetName("fLambdaPt");
 177   fLambdaPt->SetTitle("Raw #Lambda pT spectrum");
 178   fOutput->Add(fLambdaPt);
 179
 180   //----------------------
 181
 182   fLambdaBarM =
 183   new TH2F("fLambdaBarM","Mass for anti-\\Lambda", nbins, 1.065, 1.165,nbins,0.,ptMax);
 184   fLambdaBarM->GetXaxis()->SetTitle("Mass (GeV/c)");
 185   fOutput->Add(fLambdaBarM);
 186
 187   fLambdaBarSi =
 188   new TH2F("fLambdaBarSi","L_{T} vs p_{T} for anti-\\Lambda, side-band subtructed",
 189   nbins,0.,ptMax,nbins,0.,ltMax);
 190   fLambdaBarSi->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 191   fLambdaBarSi->GetYaxis()->SetTitle("L_{T} (cm)");
 192   fOutput->Add(fLambdaBarSi);
 193
 194   fLambdaBarMC =
 195   new TH2F("fLambdaBarMC","L_{T} vs p_{T} for anti-\\Lambda, from MC stack",
 196   nbins,0.,ptMax,nbins,0.,ltMax);
 197   fLambdaBarMC->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 198   fLambdaBarMC->GetYaxis()->SetTitle("L_{T} (cm)");
 199   fOutput->Add(fLambdaBarMC);
 200
 201   fLambdaBarAs =
 202   new TH2F("fLambdaBarAs","L_{T} vs p_{T} for anti-\\Lambda, associated",
 203   nbins,0.,ptMax,nbins,0.,ltMax);
 204   fLambdaBarAs->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 205   fLambdaBarAs->GetYaxis()->SetTitle("L_{T} (cm)");
 206   fOutput->Add(fLambdaBarAs);
 207
 208
 209   //----------------------
 210
 211   fLambdaBarEff=fLambdaBarAs->ProjectionX();
 212   fLambdaBarEff->SetName("fLambdaBarEff");
 213   fLambdaBarEff->SetTitle("Efficiency for anti-#Lambda");
 214   fOutput->Add(fLambdaBarEff);
 215
 216   fLambdaBarPt=fLambdaBarAs->ProjectionX();
 217   fLambdaBarPt->SetName("fLambdaBarPt");
 218   fLambdaBarPt->SetTitle("Raw anti-#Lambda pT spectrum");
 219   fOutput->Add(fLambdaBarPt);
 220
 221   //----------------------
 222
 223   fLambdaFromXi=new TH3F("fLambdaFromXi","L_{T} vs p_{T} vs p_{T} of \\Xi for \\Lambda from \\Xi",
 224   nbins,0.,ptMax,nbins,0.,ltMax,33,0.,ptMax+2);
 225   fOutput->Add(fLambdaFromXi);
 226
 227   fXiM  =
 228   new TH2F("fXiM", "\\Xi mass distribution", 50, 1.271, 1.371,33,0.,ptMax+2);
 229   fOutput->Add(fXiM);
 230
 231   fXiSiP  = new TH1F("fXiSiP", "Pt for \\Xi, side-band subracted",
 232   33,0.,ptMax+2);
 233   fOutput->Add(fXiSiP);
 234
 235
 236   fLambdaBarFromXiBar=new TH3F("fLambdaBarFromXiBar","L_{T} vs p_{T} vs p_{T} of anti-\\Xi for anti-\\Lambda from anti-\\Xi",
 237   nbins,0.,ptMax,nbins,0.,ltMax,33,0.,ptMax+2);
 238   fOutput->Add(fLambdaBarFromXiBar);
 239
 240   fXiBarM  =
 241   new TH2F("fXiBarM", "anti-\\Xi mass distribution", 50, 1.271, 1.371,33,0.,ptMax+2);
 242   fOutput->Add(fXiBarM);
 243
 244   fXiBarSiP  = new TH1F("fXiBarSiP", "Pt for anti-\\Xi, side-band subracted",
 245   33,0.,ptMax+2);
 246   fOutput->Add(fXiBarSiP);
 247
 248
 249   PostData(1, fOutput);
 250 }
 251
 252 static Bool_t AcceptTrack(const AliAODTrack *t) {
 253   if (!t->IsOn(AliAODTrack::kTPCrefit)) return kFALSE;
 254   //if (t->GetKinkIndex(0)>0) return kFALSE;
 255
 256   Float_t nCrossedRowsTPC = t->GetTPCClusterInfo(2,1);
 257   if (nCrossedRowsTPC < 70) return kFALSE;
 258   Int_t findable=t->GetTPCNclsF();
 259   if (findable <= 0) return kFALSE;
 260   if (nCrossedRowsTPC/findable < 0.8) return kFALSE;
 261
 262   if (TMath::Abs(t->Eta()) > 0.8) return kFALSE;
 263
 264   return kTRUE;
 265 }
 266
 267 Bool_t
 268 AliAnalysisTaskCTauPbPbaod::AcceptV0(const AliAODv0 *v0,const AliAODEvent *aod)
 269 {
 270   if (v0->GetOnFlyStatus()) return kFALSE;
 271
 272   Double_t cpa=v0->CosPointingAngle(aod->GetPrimaryVertex());
 273   if (cpa < fCPA) return kFALSE;
 274
 275   Double_t dca=v0->DcaV0Daughters();
 276   if (dca > fDCA) return kFALSE;
 277
 278   const AliAODTrack *ntrack=(AliAODTrack *)v0->GetDaughter(1);
 279   if (!AcceptTrack(ntrack)) return kFALSE;
 280
 281   const AliAODTrack *ptrack=(AliAODTrack *)v0->GetDaughter(0);
 282   if (!AcceptTrack(ptrack)) return kFALSE;
 283
 284   Float_t xy=v0->DcaNegToPrimVertex();
 285   if (TMath::Abs(xy)<0.1) return kFALSE;
 286   xy=v0->DcaPosToPrimVertex();
 287   if (TMath::Abs(xy)<0.1) return kFALSE;
 288
 289   Double_t xyz[3]; v0->GetSecondaryVtx(xyz);
 290   Double_t r2=xyz[0]*xyz[0] + xyz[1]*xyz[1];
 291   if (r2<0.9*0.9) return kFALSE;
 292   if (r2>100*100) return kFALSE;
 293
 294   return kTRUE;
 295 }
 296
 297 static Bool_t AcceptCascade(const AliAODcascade *cs, const AliAODEvent *aod) {
 298
 299   AliAODVertex *xi = cs->GetDecayVertexXi();
 300   const AliAODTrack *bach=(AliAODTrack *)xi->GetDaughter(0);
 301   if (!AcceptTrack(bach)) return kFALSE;
 302
 303   Double_t xy=cs->DcaBachToPrimVertex();
 304   if (TMath::Abs(xy) < 0.03) return kFALSE;
 305
 306   const AliAODVertex *vtx=aod->GetPrimaryVertex();
 307   Double_t xv=vtx->GetX(), yv=vtx->GetY(), zv=vtx->GetZ();
 308   Double_t cpa=cs->CosPointingAngleXi(xv,yv,zv);
 309   if (cpa<0.999) return kFALSE;
 310
 311   if (cs->DcaXiDaughters() > 0.3) return kFALSE;
 312
 313   return kTRUE;
 314 }
 315
 316 static Bool_t AcceptPID(const AliPIDResponse *pidResponse,
 317                         const AliAODTrack *ptrack, const TClonesArray *stack) {
 318
 319   const AliAODPid *pid=ptrack->GetDetPid();
 320   if (!pid) return kTRUE;
 321   if (pid->GetTPCmomentum() > 1.) return kTRUE;
 322
 323   if (stack) {
 324     // MC PID
 325     Int_t plab=TMath::Abs(ptrack->GetLabel());
 326     AliAODMCParticle *pp=(AliAODMCParticle*)(*stack)[plab];
 327     if (!pp) return kTRUE;
 328     if (pp->Charge() > 0) {
 329        if (pp->GetPdgCode() == kProton)    return kTRUE;
 330     } else {
 331        if (pp->GetPdgCode() == kProtonBar) return kTRUE;
 332     }
 333   } else {
 334     // Real PID
 335     Double_t nsig=pidResponse->NumberOfSigmasTPC(ptrack,AliPID::kProton);
 336     if (TMath::Abs(nsig) < 3.) return kTRUE;
 337   }
 338
 339   return kFALSE;
 340 }
 341
 342 static AliAODMCParticle*
 343 AssociateV0(const AliAODTrack *ptrack, const AliAODTrack *ntrack,
 344 const TClonesArray *stack, AliAODMCParticle *&mcp) {
 345   //
 346   // Try to associate the V0 with the daughters ptrack and ntrack
 347   // with the Monte Carlo
 348   //
 349   if (!stack) return 0;
 350
 351   Int_t nlab=TMath::Abs(ntrack->GetLabel());
 352   AliAODMCParticle *n=(AliAODMCParticle*)(*stack)[nlab];
 353   if (!n) return 0;
 354
 355   Int_t plab=TMath::Abs(ptrack->GetLabel());
 356   AliAODMCParticle *p=(AliAODMCParticle*)(*stack)[plab];
 357   if (!p) return 0;
 358
 359   Int_t imp=p->GetMother(); //V0 particle, the mother of the pos. track
 360   AliAODMCParticle *p0=(AliAODMCParticle*)(*stack)[imp];
 361   if (!p0) return 0;
 362
 363   Int_t imn=n->GetMother();
 364   if (imp != imn) {  // Check decays of the daughters
 365      return 0; // Fixme
 366   }
 367
 368   Int_t code=p0->GetPdgCode();
 369   if (code != kK0Short)
 370      if (code != kLambda0)
 371         if (code != kLambda0Bar) return 0;
 372
 373   mcp=p;
 374
 375   return p0;
 376 }
 377
 378
 379 void AliAnalysisTaskCTauPbPbaod::UserExec(Option_t *)
 380 {
 381   const Double_t yMax=0.5;
 382   const Double_t pMin=0.0;
 383   const Double_t lMax=0.001;
 384
 385   AliAODEvent *aod=(AliAODEvent *)InputEvent();
 386
 387   if (!aod) {
 388     Printf("ERROR: aod not available");
 389     return;
 390   }
 391
 392   // Vertex selection
 393   const AliAODVertex *vtx=aod->GetPrimaryVertex();
 394   if (vtx->GetNContributors()<3) return;
 395   Double_t xv=vtx->GetX(), yv=vtx->GetY(), zv=vtx->GetZ();
 396
 397   if (TMath::Abs(zv) > 10.) return ;
 398
 399
 400   // Physics selection
 401   AliAnalysisManager *mgr= AliAnalysisManager::GetAnalysisManager();
 402   AliInputEventHandler *hdr=(AliInputEventHandler*)mgr->GetInputEventHandler();
 403   UInt_t maskIsSelected = hdr->IsEventSelected();
 404   Bool_t isSelected = (maskIsSelected & AliVEvent::kMB);
 405   if (!isSelected) return;
 406
 407   fMult->Fill(-100); //event counter
 408
 409   // Centrality selection
 410   AliCentrality *cent=aod->GetCentrality();
 411   if (!cent->IsEventInCentralityClass(fCMin,fCMax,"V0M")) return;
 412
 413   AliPIDResponse *pidResponse = hdr->GetPIDResponse();
 414
 415
 416
 417   //+++++++ MC
 418   TClonesArray *stack = 0x0;
 419   Double_t mcXv=0., mcYv=0., mcZv=0.;
 420
 421   if (fIsMC) {
 422      TList *lst = aod->GetList();
 423      stack = (TClonesArray*)lst->FindObject(AliAODMCParticle::StdBranchName());
 424      if (!stack) {
 425         Printf("ERROR: stack not available");
 426         return;
 427      }
 428      AliAODMCHeader *
 429      mcHdr=(AliAODMCHeader*)lst->FindObject(AliAODMCHeader::StdBranchName());
 430
 431      mcXv=mcHdr->GetVtxX(); mcYv=mcHdr->GetVtxY(); mcZv=mcHdr->GetVtxZ();
 432
 433      Int_t ntrk=stack->GetEntriesFast();
 434      while (ntrk--) {
 435        AliAODMCParticle *p0=(AliAODMCParticle*)stack->UncheckedAt(ntrk);
 436        Int_t code=p0->GetPdgCode();
 437        if (code != kK0Short)
 438          if (code != kLambda0)
 439             if (code != kLambda0Bar) continue;
 440
 441        Int_t plab=p0->GetDaughter(0), nlab=p0->GetDaughter(1);
 442        if (nlab==plab) continue;
 443        AliAODMCParticle *part=(AliAODMCParticle*)(*stack)[plab];
 444        if (!part) continue;
 445        Double_t charge=part->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->Xv(), y=p0->Yv(), z=p0->Zv();
 453        Double_t dxmc=mcXv-x, dymc=mcYv-y, dzmc=mcZv-z;
 454        Double_t l=TMath::Sqrt(dxmc*dxmc + dymc*dymc + dzmc*dzmc);
 455
 456        if (l > lMax) continue; // secondary V0
 457
 458        x=part->Xv(); y=part->Yv();
 459        dxmc=mcXv-x; dymc=mcYv-y;
 460        Double_t lt=TMath::Sqrt(dxmc*dxmc + dymc*dymc);
 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=aod->GetNumberOfTracks();
 480   Int_t mult=0;
 481   for (Int_t i=0; i<ntrk1; i++) {
 482     AliAODTrack *t=aod->GetTrack(i);
 483     if (t->IsMuonTrack()) continue;
 484     if (!t->IsOn(AliAODTrack::kTPCrefit)) continue;
 485
 486     Double_t xyz[3];
 487     if (t->GetPosition(xyz)) continue;
 488     if (TMath::Abs(xyz[0])>3.) continue;
 489     if (TMath::Abs(xyz[1])>3.) continue;
 490
 491     Double_t pt=t->Pt(),pz=t->Pz();
 492     if (TMath::Abs(pz/pt)>0.8) continue;
 493
 494     mult++;
 495
 496     const AliAODPid *pid=t->GetDetPid();
 497     if (!pid) continue;
 498
 499     Double_t p=pid->GetTPCmomentum();
 500     Double_t dedx=pid->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 = aod->GetNumberOfV0s();
 511   while (nv0--) {
 512       AliAODv0 *v0=aod->GetV0(nv0);
 513
 514       Double_t pt=TMath::Sqrt(v0->Pt2V0());
 515       if (pt < pMin) continue;
 516
 517       if (!AcceptV0(v0,aod)) continue;
 518
 519       Double_t cpa=v0->CosPointingAngle(aod->GetPrimaryVertex());
 520       fCosPA->Fill(cpa);
 521
 522       Double_t dca=v0->DcaV0Daughters();
 523       fDtrDCA->Fill(dca);
 524
 525       const AliAODTrack *ntrack=(AliAODTrack *)v0->GetDaughter(1);
 526       const AliAODTrack *ptrack=(AliAODTrack *)v0->GetDaughter(0);
 527
 528       Double_t xyz[3]; v0->GetSecondaryVtx(xyz);
 529       Double_t dxv=xyz[0]-xv, dyv=xyz[1]-yv;
 530       Double_t lt=TMath::Sqrt(dxv*dxv + dyv*dyv);
 531
 532       if (lt/pt > 3*7.89/1.1157) continue;
 533
 534       //--- V0 switches
 535       Bool_t isK0s=kTRUE;
 536       Bool_t isLambda=kTRUE;
 537       Bool_t isLambdaBar=kTRUE;
 538
 539       if (0.4977*lt/pt > 3*2.68) isK0s=kFALSE;
 540       if (1.1157*lt/pt > 3*7.89) isLambdaBar=isLambda=kFALSE;
 541
 542       if (!AcceptPID(pidResponse, ptrack, stack)) isLambda=kFALSE;
 543       if (!AcceptPID(pidResponse, ntrack, stack)) isLambdaBar=kFALSE;
 544
 545       Double_t yK0s=TMath::Abs(v0->RapK0Short());
 546       Double_t yLam=TMath::Abs(v0->RapLambda());
 547       if (yK0s > yMax) isK0s=kFALSE;
 548       if (yLam > yMax) isLambda=isLambdaBar=kFALSE;
 549       //---
 550
 551       Double_t mass=0., m=0., s=0.;
 552       if (isK0s) {
 553          mass=v0->MassK0Short();
 554          fK0sM->Fill(mass,pt);
 555
 556          m=TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass();
 557          s=0.0044 + (0.008-0.0044)/(10-1)*(pt - 1.);
 558          if (TMath::Abs(m-mass) < 3*s) {
 559             fK0sSi->Fill(pt,lt);
 560          } else {
 561             isK0s=kFALSE;
 562          }
 563          if (TMath::Abs(m-mass + 4.5*s) < 1.5*s) {
 564             fK0sSi->Fill(pt,lt,-1);
 565          }
 566          if (TMath::Abs(m-mass - 4.5*s) < 1.5*s) {
 567             fK0sSi->Fill(pt,lt,-1);
 568          }
 569       }
 570
 571       if (isLambda) {
 572          mass=v0->MassLambda();
 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          mass=v0->MassAntiLambda();
 594          fLambdaBarM->Fill(mass,pt);
 595
 596          m=TDatabasePDG::Instance()->GetParticle(kLambda0Bar)->Mass();
 597          //s=0.0027 + (0.004-0.0027)/(10-1)*(pt-1);
 598          //s=0.0015 + (0.002-0.0015)/(2.6-1)*(pt-1);
 599          s=0.0023 + (0.004-0.0023)/(6-1)*(pt-1);
 600          if (TMath::Abs(m-mass) < 3*s) {
 601             fLambdaBarSi->Fill(pt,lt);
 602          } else {
 603             isLambdaBar=kFALSE;
 604          }
 605          if (TMath::Abs(m-mass + 4.5*s) < 1.5*s) {
 606             fLambdaBarSi->Fill(pt,lt,-1);
 607          }
 608          if (TMath::Abs(m-mass - 4.5*s) < 1.5*s) {
 609             fLambdaBarSi->Fill(pt,lt,-1);
 610          }
 611       }
 612
 613       if (!fIsMC) continue;
 614
 615       //++++++ MC
 616       if (!isK0s)
 617          if (!isLambda)
 618              if (!isLambdaBar) continue;//check MC only for the accepted V0s
 619
 620       AliAODMCParticle *mcp=0;
 621       AliAODMCParticle *mc0=AssociateV0(ptrack,ntrack,stack,mcp);
 622       if (!mc0) continue;
 623
 624       Double_t ptAs=mc0->Pt();
 625       if (ptAs < pMin) continue;
 626       Double_t yAs=mc0->Y();
 627       if (TMath::Abs(yAs) > yMax) continue;
 628
 629       Int_t code=mc0->GetPdgCode();
 630
 631       Double_t dx = mcXv - mcp->Xv(), dy = mcYv - mcp->Yv();
 632       Double_t ltAs=TMath::Sqrt(dx*dx + dy*dy);
 633
 634       Double_t dz=mcZv - mc0->Zv(); dy=mcYv - mc0->Yv(); dx=mcXv - mc0->Xv();
 635       Double_t l = TMath::Sqrt(dx*dx + dy*dy + dz*dz);
 636       if (l<lMax) { // Primary V0s
 637          switch (code) {
 638          case kK0Short:
 639             if (isK0s)       fK0sAs->Fill(ptAs,ltAs);
 640             break;
 641          case kLambda0:
 642             if (isLambda)    fLambdaAs->Fill(ptAs,ltAs);
 643             break;
 644          case kLambda0Bar:
 645             if (isLambdaBar) fLambdaBarAs->Fill(ptAs,ltAs);
 646             break;
 647          default: break;
 648          }
 649       } else {
 650          if (code==kK0Short) continue;
 651
 652          Int_t nx=mc0->GetMother();
 653          AliAODMCParticle *xi=(AliAODMCParticle*)(*stack)[nx];
 654          if (!xi) continue;
 655          Int_t xcode=xi->GetPdgCode();
 656
 657          switch (code) {
 658          case kLambda0:
 659             if (isLambda)
 660             if ((xcode==kXiMinus) || (xcode==3322))
 661                fLambdaFromXi->Fill(ptAs,ltAs,xi->Pt());
 662             break;
 663          case kLambda0Bar:
 664             if (isLambdaBar)
 665             if ((xcode==kXiPlusBar)||(xcode==-3322))
 666                fLambdaBarFromXiBar->Fill(ptAs,ltAs,xi->Pt());
 667             break;
 668          default: break;
 669          }
 670       }
 671       //++++++
 672
 673   }
 674
 675   Int_t ncs=aod->GetNumberOfCascades();
 676   for (Int_t i=0; i<ncs; i++) {
 677       AliAODcascade *cs=aod->GetCascade(i);
 678
 679       Double_t pt=TMath::Sqrt(cs->Pt2Xi());
 680       if (pt < pMin) continue;
 681       if (TMath::Abs(cs->RapXi()) > yMax) continue;
 682       if (!AcceptCascade(cs,aod)) continue;
 683
 684       const AliAODv0 *v0=(AliAODv0*)cs;
 685       if (v0->RapLambda() > yMax) continue;
 686       if (!AcceptV0(v0,aod)) continue;
 687
 688       //--- Cascade switches
 689       Bool_t isXiMinus=kTRUE;
 690       Bool_t isXiPlusBar=kTRUE;
 691
 692       const AliAODTrack *ptrack=(AliAODTrack *)v0->GetDaughter(0);
 693       if (!AcceptPID(pidResponse, ptrack, stack)) isXiMinus=kFALSE;
 694
 695       const AliAODTrack *ntrack=(AliAODTrack *)v0->GetDaughter(1);
 696       if (!AcceptPID(pidResponse, ntrack, stack)) isXiPlusBar=kFALSE;
 697
 698       Int_t charge=cs->ChargeXi();
 699       if (charge > 0) isXiMinus=kFALSE;
 700       if (charge < 0) isXiPlusBar=kFALSE;
 701       //---
 702
 703       if (isXiMinus) {
 704          Double_t mass=cs->MassXi();
 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          Double_t mass=cs->MassXi();
 724          fXiBarM->Fill(mass,pt);
 725          Double_t m=TDatabasePDG::Instance()->GetParticle(kXiPlusBar)->Mass();
 726          //Double_t s=0.0037;
 727          Double_t s=0.002 + (0.0032-0.002)/(6-1.5)*(pt-1.5);
 728          if (TMath::Abs(m-mass) < 3*s) {
 729             fXiBarSiP->Fill(pt);
 730          } else {
 731             isXiPlusBar=kFALSE;
 732          }
 733          if (TMath::Abs(m-mass + 4.5*s) < 1.5*s) {
 734             fXiBarSiP->Fill(pt,-1);
 735          }
 736          if (TMath::Abs(m-mass - 4.5*s) < 1.5*s) {
 737             fXiBarSiP->Fill(pt,-1);
 738          }
 739       }
 740
 741       if (!fIsMC) continue;
 742
 743       //++++++ MC
 744       if (!isXiMinus)
 745          if (!isXiPlusBar) continue;//check MC only for the accepted cascades
 746       // Here is the future association with MC
 747
 748   }
 749
 750 }
 751
 752 void AliAnalysisTaskCTauPbPbaod::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 }