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