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