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 static Int_t    nbins=100;  // number of bins
  35 static Double_t lMin=0.0, lMax=100.;
  36 static Double_t pMin=0.0, pMax=10.;
  37 static Double_t yMax=0.5;
  38
  39
  40 //
  41 //  This is a little task for checking the c*tau of the strange particles
  42 //
  43
  44 AliAnalysisTaskCTauPbPbaod::AliAnalysisTaskCTauPbPbaod(const char *name) :
  45 AliAnalysisTaskSE(name),
  46 fIsMC(kFALSE),
  47 fCMin(0.),
  48 fCMax(90.),
  49 fOutput(0),
  50 fMult(0),
  51 fdEdx(0),
  52 fdEdxPid(0),
  53
  54 fK0sM(0),
  55 fK0sSi(0),
  56 fK0sMC(0),
  57 fK0sAs(0),
  58
  59 fLambdaM(0),
  60 fLambdaSi(0),
  61 fLambdaMC(0),
  62 fLambdaAs(0),
  63
  64 fCPA(0),
  65 fDCA(0),
  66
  67 fLambdaEff(0),
  68 fLambdaPt(0),
  69
  70 fLambdaFromXi(0),
  71 fXiM(0),
  72 fXiSiP(0)
  73 {
  74   // Constructor. Initialization of pointers
  75   DefineOutput(1, TList::Class());
  76 }
  77
  78 void AliAnalysisTaskCTauPbPbaod::UserCreateOutputObjects()
  79 {
  80   fOutput = new TList();
  81   fOutput->SetOwner();
  82
  83
  84   fMult=new TH1F("fMult","Multiplicity",1100,0.,3300);
  85   fMult->GetXaxis()->SetTitle("N tracks");
  86   fOutput->Add(fMult);
  87
  88   fdEdx=new TH2F("fdEdx","dE/dx",50,0.2,3,50,0.,6.);
  89   fOutput->Add(fdEdx);
  90
  91   fdEdxPid=new TH2F("fdEdxPid","dE/dx with PID",50,0.2,3,50,0.,6.);
  92   fOutput->Add(fdEdxPid);
  93
  94   fK0sM =
  95   new TH2F("fK0sM", "Mass for K^{0}_{s}", nbins/2,0.448,0.548,nbins,pMin,pMax);
  96   fK0sM->GetXaxis()->SetTitle("Mass (GeV/c)");
  97   fOutput->Add(fK0sM);
  98
  99   fK0sSi =
 100   new TH2F("fK0sSi","L_{T} vs p_{T} for K^{0}_{s}, side-band subtracted",
 101   nbins,pMin,pMax,nbins,lMin,lMax);
 102   fK0sSi->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 103   fK0sSi->GetYaxis()->SetTitle("L_{T} (cm)");
 104   fOutput->Add(fK0sSi);
 105
 106   fK0sMC =
 107   new TH2F("fK0sMC","L_{T} vs p_{T} for K^{0}_{s}, from MC stack",
 108   nbins,pMin,pMax,nbins,lMin,lMax);
 109   fK0sMC->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 110   fK0sMC->GetYaxis()->SetTitle("L_{T} (cm)");
 111   fOutput->Add(fK0sMC);
 112
 113   fK0sAs =
 114   new TH2F("fK0sAs", "L_{T} vs p_{T} for K^{0}_{s}, associated",
 115   nbins,pMin,pMax,nbins,lMin,lMax);
 116   fK0sAs->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 117   fK0sAs->GetYaxis()->SetTitle("L_{T} (cm)");
 118   fOutput->Add(fK0sAs);
 119
 120   //----------------------
 121
 122   fLambdaM =
 123   new TH2F("fLambdaM","Mass for \\Lambda", nbins, 1.065, 1.165,nbins,pMin,pMax);
 124   fLambdaM->GetXaxis()->SetTitle("Mass (GeV/c)");
 125   fOutput->Add(fLambdaM);
 126
 127   fLambdaSi =
 128   new TH2F("fLambdaSi","L_{T} vs p_{T} for \\Lambda, side-band subtructed",
 129   nbins,pMin,pMax,nbins,lMin,lMax);
 130   fLambdaSi->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 131   fLambdaSi->GetYaxis()->SetTitle("L_{T} (cm)");
 132   fOutput->Add(fLambdaSi);
 133
 134   fLambdaMC =
 135   new TH2F("fLambdaMC","L_{T} vs p_{T} for \\Lambda, from MC stack",
 136   nbins,pMin,pMax,nbins,lMin,lMax);
 137   fLambdaMC->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 138   fLambdaMC->GetYaxis()->SetTitle("L_{T} (cm)");
 139   fOutput->Add(fLambdaMC);
 140
 141   fLambdaAs =
 142   new TH2F("fLambdaAs","L_{T} vs p_{T} for \\Lambda, associated",
 143   nbins,pMin,pMax,nbins,lMin,lMax);
 144   fLambdaAs->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 145   fLambdaAs->GetYaxis()->SetTitle("L_{T} (cm)");
 146   fOutput->Add(fLambdaAs);
 147
 148   fCPA=new TH1F("fCPA","Cosine of the pointing angle",30,0.9978,1.);
 149   fOutput->Add(fCPA);
 150   fDCA=new TH1F("fDCA","DCA between the daughters",30,0.,1.1);
 151   fOutput->Add(fDCA);
 152
 153   fLambdaEff=fLambdaAs->ProjectionX();
 154   fLambdaEff->SetName("fLambdaEff");
 155   fLambdaEff->SetTitle("Efficiency for #Lambda");
 156   fOutput->Add(fLambdaEff);
 157
 158   fLambdaPt=fLambdaAs->ProjectionX();
 159   fLambdaPt->SetName("fLambdaPt");
 160   fLambdaPt->SetTitle("Raw #Lambda pT spectrum");
 161   fOutput->Add(fLambdaPt);
 162
 163   //----------------------
 164
 165   fLambdaFromXi=new TH3F("fLambdaFromXi","L_{T} vs p_{T} vs p_{T} of \\Xi for \\Lambda from Xi",
 166   nbins,pMin,pMax,nbins,lMin,lMax,33,pMin,pMax+2);
 167   fOutput->Add(fLambdaFromXi);
 168
 169   fXiM  =
 170   new TH2F("fXiM", "\\Xi mass distribution", 50, 1.271, 1.371,12,pMin,pMax+2);
 171   fOutput->Add(fXiM);
 172
 173   fXiSiP  = new TH1F("fXiSiP", "Pt for \\Xi, side-band subracted",
 174   33,pMin,pMax+2);
 175   fOutput->Add(fXiSiP);
 176
 177
 178   PostData(1, fOutput);
 179 }
 180
 181 static Bool_t AcceptTrack(const AliAODTrack *t) {
 182   if (!t->IsOn(AliAODTrack::kTPCrefit)) return kFALSE;
 183   //if (t->GetKinkIndex(0)>0) return kFALSE;
 184
 185   Float_t nCrossedRowsTPC = t->GetTPCClusterInfo(2,1);
 186   if (nCrossedRowsTPC < 70) return kFALSE;
 187   Int_t findable=t->GetTPCNclsF();
 188   if (findable <= 0) return kFALSE;
 189   if (nCrossedRowsTPC/findable < 0.8) return kFALSE;
 190
 191   return kTRUE;
 192 }
 193
 194 static Bool_t AcceptV0(const AliAODv0 *v1, const AliAODEvent *aod) {
 195
 196   if (v1->GetOnFlyStatus()) return kFALSE;
 197
 198   if (v1->Pt() < pMin) return kFALSE;
 199
 200   const AliAODTrack *ntrack1=(AliAODTrack *)v1->GetDaughter(1);
 201   if (!AcceptTrack(ntrack1)) return kFALSE;
 202
 203   const AliAODTrack *ptrack1=(AliAODTrack *)v1->GetDaughter(0);
 204   if (!AcceptTrack(ptrack1)) return kFALSE;
 205
 206   Float_t xy=v1->DcaNegToPrimVertex();
 207   if (TMath::Abs(xy)<0.1) return kFALSE;
 208   xy=v1->DcaPosToPrimVertex();
 209   if (TMath::Abs(xy)<0.1) return kFALSE;
 210
 211   Double_t dca=v1->DcaV0Daughters();
 212   if (dca>1.0) return kFALSE;
 213   //if (dca>0.7) return kFALSE;
 214   //if (dca>0.4) return kFALSE;
 215
 216   Double_t cpa=v1->CosPointingAngle(aod->GetPrimaryVertex());
 217   if (cpa<0.998) return kFALSE;
 218   //if (cpa<0.99875) return kFALSE;
 219   //if (cpa<0.9995) return kFALSE;
 220
 221   Double_t xyz[3]; v1->GetSecondaryVtx(xyz);
 222   Double_t r2=xyz[0]*xyz[0] + xyz[1]*xyz[1];
 223   if (r2<0.9*0.9) return kFALSE;
 224   if (r2>100*100) return kFALSE;
 225
 226   return kTRUE;
 227 }
 228
 229 static Bool_t AcceptCascade(const AliAODcascade *cs, const AliAODEvent *aod) {
 230
 231   if (cs->Pt() < pMin) return kFALSE;
 232
 233   AliAODVertex *xi = cs->GetDecayVertexXi();
 234   const AliAODTrack *bach=(AliAODTrack *)xi->GetDaughter(0);
 235   if (!AcceptTrack(bach)) return kFALSE;
 236
 237   Double_t xy=cs->DcaBachToPrimVertex();
 238   if (TMath::Abs(xy) < 0.03) return kFALSE;
 239
 240   const AliAODVertex *vtx=aod->GetPrimaryVertex();
 241   Double_t xv=vtx->GetX(), yv=vtx->GetY(), zv=vtx->GetZ();
 242   Double_t cpa=cs->CosPointingAngleXi(xv,yv,zv);
 243   if (cpa<0.999) return kFALSE;
 244
 245   if (cs->DcaXiDaughters() > 0.3) return kFALSE;
 246
 247   return kTRUE;
 248 }
 249
 250 void AliAnalysisTaskCTauPbPbaod::UserExec(Option_t *)
 251 {
 252
 253   AliAODEvent *aod=(AliAODEvent *)InputEvent();
 254
 255   if (!aod) {
 256     Printf("ERROR: aod not available");
 257     return;
 258   }
 259
 260   fMult->Fill(-100); //event counter
 261
 262   // Physics selection
 263   AliAnalysisManager *mgr= AliAnalysisManager::GetAnalysisManager();
 264   AliInputEventHandler *hdr=(AliInputEventHandler*)mgr->GetInputEventHandler();
 265   UInt_t maskIsSelected = hdr->IsEventSelected();
 266   Bool_t isSelected = (maskIsSelected & AliVEvent::kMB);
 267   if (!isSelected) return;
 268
 269   // Centrality selection
 270   AliCentrality *cent=aod->GetCentrality();
 271   if (!cent->IsEventInCentralityClass(fCMin,fCMax,"V0M")) return;
 272
 273   const AliAODVertex *vtx=aod->GetPrimaryVertex();
 274   if (vtx->GetNContributors()<3) return;
 275
 276   Double_t xv=vtx->GetX(), yv=vtx->GetY(), zv=vtx->GetZ();
 277
 278   if (TMath::Abs(zv) > 10.) return ;
 279
 280   AliPIDResponse *pidResponse = hdr->GetPIDResponse();
 281
 282   //fMult->Fill(-100); //event counter
 283
 284   //+++++++ MC
 285   TClonesArray *stack = 0x0;
 286   Double_t mcXv=0., mcYv=0., mcZv=0.;
 287
 288   if (fIsMC) {
 289      TList *lst = aod->GetList();
 290      stack = (TClonesArray*)lst->FindObject(AliAODMCParticle::StdBranchName());
 291      if (!stack) {
 292         Printf("ERROR: stack not available");
 293         return;
 294      }
 295      AliAODMCHeader *
 296      mcHdr=(AliAODMCHeader*)lst->FindObject(AliAODMCHeader::StdBranchName());
 297
 298      mcXv=mcHdr->GetVtxX(); mcYv=mcHdr->GetVtxY(); mcZv=mcHdr->GetVtxZ();
 299
 300      Int_t ntrk=stack->GetEntriesFast(), ntrk0=ntrk;
 301      while (ntrk--) {
 302        AliAODMCParticle *p0=(AliAODMCParticle*)stack->UncheckedAt(ntrk);
 303        Int_t code=p0->GetPdgCode();
 304        if (code != kK0Short)
 305          if (code != kLambda0) continue;
 306
 307        Int_t plab=p0->GetDaughter(0), nlab=p0->GetDaughter(1);
 308        if (nlab==plab) continue;
 309        if (nlab<0) continue;
 310        if (plab<0) continue;
 311        if (nlab>=ntrk0) continue;
 312        if (plab>=ntrk0) continue;
 313        AliAODMCParticle *part = (AliAODMCParticle *)stack->UncheckedAt(plab);
 314        if (!part) continue;
 315        Double_t charge=part->Charge();
 316        if (charge == 0.) continue;
 317
 318        Double_t pt=p0->Pt();
 319        if (pt<pMin) continue;
 320        if (TMath::Abs(p0->Y())>yMax) continue;
 321
 322        Double_t x=p0->Xv(), y=p0->Yv(), z=p0->Zv();
 323        Double_t dx=mcXv-x, dy=mcYv-y, dz=mcZv-z;
 324        Double_t l=TMath::Sqrt(dx*dx + dy*dy + dz*dz);
 325
 326        if (l > 0.01) continue; // secondary V0
 327
 328        x=part->Xv(); y=part->Yv();
 329        dx=mcXv-x; dy=mcYv-y;
 330        Double_t lt=TMath::Sqrt(dx*dx + dy*dy);
 331
 332        if (code == kK0Short) {
 333           fK0sMC->Fill(pt,lt);
 334        }
 335        if (code == kLambda0) {
 336           fLambdaMC->Fill(pt,lt);
 337        }
 338      }
 339   }
 340
 341
 342   Int_t ntrk=aod->GetNumberOfTracks();
 343   Int_t mult=0;
 344   Double_t nsig;
 345   for (Int_t i=0; i<ntrk; i++) {
 346     AliAODTrack *t=aod->GetTrack(i);
 347     if (t->IsMuonTrack()) continue;
 348     if (!t->IsOn(AliAODTrack::kTPCrefit)) continue;
 349
 350     Double_t xyz[3];
 351     if (t->GetPosition(xyz)) continue;
 352     if (TMath::Abs(xyz[0])>3.) continue;
 353     if (TMath::Abs(xyz[1])>3.) continue;
 354
 355     Double_t pt=t->Pt(),pz=t->Pz();
 356     if (TMath::Abs(pz/pt)>0.8) continue;
 357
 358     mult++;
 359
 360     const AliAODPid *pid=t->GetDetPid();
 361     if (!pid) continue;
 362
 363     Double_t p=pid->GetTPCmomentum();
 364     Double_t dedx=pid->GetTPCsignal()/47.;
 365     fdEdx->Fill(p,dedx,1);
 366
 367     nsig=pidResponse->NumberOfSigmasTPC(t,AliPID::kProton);
 368     if (TMath::Abs(nsig) < 3.) fdEdxPid->Fill(p,dedx,1);
 369
 370   }
 371   fMult->Fill(mult);
 372
 373
 374   Int_t nv0 = aod->GetNumberOfV0s();
 375   while (nv0--) {
 376       AliAODv0 *v0=aod->GetV0(nv0);
 377
 378       if (!AcceptV0(v0,aod)) continue;
 379
 380       const AliAODTrack *ntrack=(AliAODTrack *)v0->GetDaughter(1);
 381       const AliAODTrack *ptrack=(AliAODTrack *)v0->GetDaughter(0);
 382
 383       Double_t xyz[3]; v0->GetSecondaryVtx(xyz);
 384       Double_t dx=xyz[0]-xv, dy=xyz[1]-yv;
 385       Double_t lt=TMath::Sqrt(dx*dx + dy*dy);
 386
 387       Double_t pt=TMath::Sqrt(v0->Pt2V0());
 388
 389       Bool_t ctK=kTRUE; if (0.4977*lt/pt > 3*2.68) ctK=kFALSE;
 390       Bool_t ctL=kTRUE; if (1.1157*lt/pt > 3*7.89) ctL=kFALSE;
 391
 392       //+++++++ MC
 393       if (stack) {
 394          Int_t ntrk=stack->GetEntriesFast();
 395
 396          Int_t nlab=TMath::Abs(ntrack->GetLabel());
 397          Int_t plab=TMath::Abs(ptrack->GetLabel());
 398
 399          if (nlab<0) goto noas;
 400          if (nlab>=ntrk) goto noas;
 401          if (plab<0) goto noas;
 402          if (plab>=ntrk) goto noas;
 403
 404          AliAODMCParticle *np=(AliAODMCParticle*)stack->UncheckedAt(nlab);
 405          AliAODMCParticle *pp=(AliAODMCParticle*)stack->UncheckedAt(plab);
 406
 407          Int_t i0=pp->GetMother();
 408          //if (np->GetMother() != i0) goto noas;
 409
 410          Int_t in0=np->GetMother();
 411          if (in0<0) goto noas;
 412          if (in0>=ntrk) goto noas;
 413          if (in0 != i0) { // did the negative daughter decay ?
 414             AliAODMCParticle *nnp=(AliAODMCParticle*)stack->UncheckedAt(in0);
 415             if (nnp->GetMother() != i0) goto noas;
 416          }
 417
 418          if (i0<0) goto noas;
 419          if (i0>=ntrk) goto noas;
 420          AliAODMCParticle *p0=(AliAODMCParticle*)stack->UncheckedAt(i0);
 421
 422          Int_t code=p0->GetPdgCode();
 423          if (code != kK0Short)
 424             if (code != kLambda0) goto noas;
 425
 426          if (p0->Pt()<pMin) goto noas;
 427          if (TMath::Abs(p0->Y())>yMax ) goto noas;
 428
 429          Double_t dz=mcZv - p0->Zv(), dy=mcYv - p0->Yv(), dx=mcXv - p0->Xv();
 430          Double_t l = TMath::Sqrt(dx*dx + dy*dy + dz*dz);
 431
 432          dx = mcXv - pp->Xv(); dy = mcYv - pp->Yv();
 433          Double_t ltAs=TMath::Sqrt(dx*dx + dy*dy);
 434          Double_t ptAs=p0->Pt();
 435
 436          if (l > 0.01) { // Secondary V0
 437            if (code != kLambda0) goto noas;
 438            Int_t nx=p0->GetMother();
 439            if (nx<0) goto noas;
 440            if (nx>=ntrk) goto noas;
 441            AliAODMCParticle *xi=(AliAODMCParticle*)stack->UncheckedAt(nx);
 442            Int_t xcode=xi->GetPdgCode();
 443            if ( xcode != kXiMinus )
 444              if( xcode != 3322 ) goto noas;
 445            fLambdaFromXi->Fill(ptAs,ltAs,xi->Pt());
 446          } else {
 447            if (code == kLambda0) {
 448               if (ctL) fLambdaAs->Fill(ptAs,ltAs);
 449            }else {
 450               if (ctK) fK0sAs->Fill(ptAs,ltAs);
 451            }
 452          }
 453
 454       }
 455       //++++++++
 456
 457   noas:
 458
 459       Double_t dca=v0->DcaV0Daughters();
 460       Double_t cpa=v0->CosPointingAngle(aod->GetPrimaryVertex());
 461
 462       Double_t mass=0., m=0., s=0.;
 463       if (ctK)
 464       if (TMath::Abs(v0->RapK0Short())<yMax) {
 465          mass=v0->MassK0Short();
 466          fK0sM->Fill(mass,pt);
 467
 468          m=TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass();
 469          s=0.0044 + (0.008-0.0044)/(10-1)*(pt - 1.);
 470          if (TMath::Abs(m-mass) < 3*s) {
 471             fK0sSi->Fill(pt,lt);
 472          }
 473          if (TMath::Abs(m-mass + 4.5*s) < 1.5*s) {
 474             fK0sSi->Fill(pt,lt,-1);
 475          }
 476          if (TMath::Abs(m-mass - 4.5*s) < 1.5*s) {
 477             fK0sSi->Fill(pt,lt,-1);
 478          }
 479       }
 480
 481       if (ctL)
 482       if (TMath::Abs(v0->RapLambda())<yMax) {
 483          const AliAODPid *pid=ptrack->GetDetPid();
 484          if (pid) {
 485             Double_t p=pid->GetTPCmomentum();
 486             if (p<1.) {
 487                nsig=pidResponse->NumberOfSigmasTPC(ptrack,AliPID::kProton);
 488                if (TMath::Abs(nsig) > 3.) continue;
 489             }
 490          }
 491          mass=v0->MassLambda();
 492          fLambdaM->Fill(mass,pt);
 493
 494          m=TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass();
 495          //s=0.0027 + (0.004-0.0027)/(10-1)*(pt-1);
 496          //s=0.0015 + (0.002-0.0015)/(2.6-1)*(pt-1);
 497          s=0.0023 + (0.004-0.0023)/(6-1)*(pt-1);
 498          if (TMath::Abs(m-mass) < 3*s) {
 499             fLambdaSi->Fill(pt,lt);
 500             fCPA->Fill(cpa,1);
 501             fDCA->Fill(dca,1);
 502          }
 503          if (TMath::Abs(m-mass + 4.5*s) < 1.5*s) {
 504             fLambdaSi->Fill(pt,lt,-1);
 505             fCPA->Fill(cpa,-1);
 506             fDCA->Fill(dca,-1);
 507          }
 508          if (TMath::Abs(m-mass - 4.5*s) < 1.5*s) {
 509             fLambdaSi->Fill(pt,lt,-1);
 510             fCPA->Fill(cpa,-1);
 511             fDCA->Fill(dca,-1);
 512          }
 513       }
 514   }
 515
 516   Int_t ncs=aod->GetNumberOfCascades();
 517   for (Int_t i=0; i<ncs; i++) {
 518       AliAODcascade *cs=aod->GetCascade(i);
 519
 520       if (TMath::Abs(cs->RapXi()) > yMax) continue;
 521       if (!AcceptCascade(cs,aod)) continue;
 522
 523       const AliAODv0 *v0=(AliAODv0*)cs;
 524       if (v0->RapLambda() > yMax) continue;
 525       if (!AcceptV0(v0,aod)) continue;
 526
 527       Double_t pt=TMath::Sqrt(cs->Pt2Xi());
 528
 529       Int_t charge=cs->ChargeXi();
 530       if (charge < 0) {
 531          //AliAODVertex *xi = cs->GetDecayVertexXi();
 532          //const AliAODv0 *v0=dynamic_cast<AliAODv0 *>(xi->GetDaughter(1));
 533          const AliAODTrack *ptrack=(AliAODTrack *)v0->GetDaughter(0);
 534          const AliAODPid *pid=ptrack->GetDetPid();
 535          if (pid) {
 536             Double_t p=pid->GetTPCmomentum();
 537             if (p<1.) {
 538                nsig=pidResponse->NumberOfSigmasTPC(ptrack,AliPID::kProton);
 539                if (TMath::Abs(nsig) > 3.) continue;
 540             }
 541          }
 542          Double_t mass=cs->MassXi();
 543          fXiM->Fill(mass,pt);
 544          Double_t m=TDatabasePDG::Instance()->GetParticle(kXiMinus)->Mass();
 545          //Double_t s=0.0037;
 546          Double_t s=0.002 + (0.0032-0.002)/(6-1.5)*(pt-1.5);
 547          if (TMath::Abs(m-mass) < 3*s) {
 548             fXiSiP->Fill(pt);
 549          }
 550          if (TMath::Abs(m-mass + 4.5*s) < 1.5*s) {
 551             fXiSiP->Fill(pt,-1);
 552          }
 553          if (TMath::Abs(m-mass - 4.5*s) < 1.5*s) {
 554             fXiSiP->Fill(pt,-1);
 555          }
 556       }
 557   }
 558
 559 }
 560
 561 void AliAnalysisTaskCTauPbPbaod::Terminate(Option_t *)
 562 {
 563    // The Terminate() function is the last function to be called during
 564    // a query. It always runs on the client, it can be used to present
 565    // the results graphically or save the results to file.
 566
 567   fOutput=(TList*)GetOutputData(1);
 568   if (!fOutput) {
 569      Printf("ERROR: fOutput not available");
 570      return;
 571   }
 572
 573   fMult = dynamic_cast<TH1F*>(fOutput->FindObject("fMult")) ;
 574   if (!fMult) {
 575      Printf("ERROR: fMult not available");
 576      return;
 577   }
 578
 579   fdEdx = dynamic_cast<TH2F*>(fOutput->FindObject("fdEdx")) ;
 580   if (!fdEdx) {
 581      Printf("ERROR: fdEdx not available");
 582      return;
 583   }
 584
 585   fdEdxPid = dynamic_cast<TH2F*>(fOutput->FindObject("fdEdxPid")) ;
 586   if (!fdEdxPid) {
 587      Printf("ERROR: fdEdxPid not available");
 588      return;
 589   }
 590
 591
 592   fK0sMC = dynamic_cast<TH2F*>(fOutput->FindObject("fK0sMC")) ;
 593   if (!fK0sMC) {
 594      Printf("ERROR: fK0sMC not available");
 595      return;
 596   }
 597   TH1D *k0sMcPx=fK0sMC->ProjectionX(); k0sMcPx->Sumw2();
 598   fK0sAs = dynamic_cast<TH2F*>(fOutput->FindObject("fK0sAs")) ;
 599   if (!fK0sAs) {
 600      Printf("ERROR: fK0sAs not available");
 601      return;
 602   }
 603   TH1D *k0sAsPx=fK0sAs->ProjectionX();
 604   k0sAsPx->Sumw2(); //k0sAsPx->Scale(0.69);
 605
 606
 607
 608   fLambdaFromXi = dynamic_cast<TH3F*>(fOutput->FindObject("fLambdaFromXi")) ;
 609   if (!fLambdaFromXi) {
 610      Printf("ERROR: fLambdaFromXi not available");
 611      return;
 612   }
 613   TH1D *lambdaFromXiPx=fLambdaFromXi->ProjectionX(); lambdaFromXiPx->Sumw2();
 614
 615
 616   fLambdaMC = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaMC")) ;
 617   if (!fLambdaMC) {
 618      Printf("ERROR: fLambdaMC not available");
 619      return;
 620   }
 621   TH1D *lambdaMcPx=fLambdaMC->ProjectionX(); lambdaMcPx->Sumw2();
 622
 623   fLambdaAs = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaAs")) ;
 624   if (!fLambdaAs) {
 625      Printf("ERROR: fLambdaAs not available");
 626      return;
 627   }
 628   TH1D *lambdaAsPx=fLambdaAs->ProjectionX();
 629   lambdaAsPx->Sumw2(); //lambdaAsPx->Scale(0.64);
 630
 631   fLambdaSi = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaSi")) ;
 632   if (!fLambdaSi) {
 633      Printf("ERROR: fLambdaSi not available");
 634      return;
 635   }
 636   TH1D *lambdaSiPx=fLambdaSi->ProjectionX();
 637   lambdaSiPx->SetName("fLambdaPt");
 638   lambdaSiPx->Sumw2();
 639
 640   fLambdaEff = dynamic_cast<TH1D*>(fOutput->FindObject("fLambdaEff")) ;
 641   if (!fLambdaEff) {
 642      Printf("ERROR: fLambdaEff not available");
 643      return;
 644   }
 645   fLambdaPt = dynamic_cast<TH1D*>(fOutput->FindObject("fLambdaPt")) ;
 646   if (!fLambdaPt) {
 647      Printf("ERROR: fLambdaPt not available");
 648      return;
 649   }
 650
 651
 652   if (!gROOT->IsBatch()) {
 653
 654     TCanvas *c1 = new TCanvas("c1","Mulitplicity");
 655     c1->SetLogy();
 656     fMult->DrawCopy() ;
 657
 658     new TCanvas("c2","dE/dx");
 659     fdEdx->DrawCopy() ;
 660
 661     new TCanvas("c3","dE/dx with PID");
 662     fdEdxPid->DrawCopy() ;
 663
 664     if (fIsMC) {
 665        /*
 666        TH1D effK(*k0sAsPx); effK.SetTitle("Efficiency for K0s");
 667        effK.Divide(k0sAsPx,k0sMcPx,1,1,"b");
 668        new TCanvas("c4","Efficiency for K0s");
 669        effK.DrawCopy("E") ;
 670        */
 671
 672        fLambdaEff->Divide(lambdaAsPx,lambdaMcPx,1,1,"b");
 673        new TCanvas("c5","Efficiency for #Lambda");
 674        fLambdaEff->DrawCopy("E") ;
 675
 676        lambdaSiPx->Add(lambdaFromXiPx,-1);
 677        lambdaSiPx->Divide(fLambdaEff);
 678
 679        new TCanvas("c6","Corrected #Lambda pt");
 680        lambdaSiPx->SetTitle("Corrected #Lambda pt");
 681       *fLambdaPt = *lambdaSiPx;
 682        fLambdaPt->SetLineColor(2);
 683        fLambdaPt->DrawCopy("E");
 684
 685        lambdaMcPx->DrawCopy("same");
 686
 687     } else {
 688        new TCanvas("c6","Raw #Lambda pt");
 689        lambdaSiPx->SetTitle("Raw #Lambda pt");
 690       *fLambdaPt = *lambdaSiPx;
 691        fLambdaPt->SetLineColor(2);
 692        fLambdaPt->DrawCopy("E");
 693     }
 694   }
 695 }