PWGLF/STRANGENESS/LambdaK0PbPb/AliAnalysisTaskCTauPbPbaod.cxx

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