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