An example of a c*tau analysis using ESDs
authorbelikov <belikov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 4 Oct 2011 10:35:26 +0000 (10:35 +0000)
committerbelikov <belikov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 4 Oct 2011 10:35:26 +0000 (10:35 +0000)
PWG2/SPECTRA/LambdaK0PbPb/AddTaskCTau.C [new file with mode: 0644]
PWG2/SPECTRA/LambdaK0PbPb/AliAnalysisTaskCTauPbPb.cxx [new file with mode: 0644]
PWG2/SPECTRA/LambdaK0PbPb/AliAnalysisTaskCTauPbPb.h [new file with mode: 0644]

diff --git a/PWG2/SPECTRA/LambdaK0PbPb/AddTaskCTau.C b/PWG2/SPECTRA/LambdaK0PbPb/AddTaskCTau.C
new file mode 100644 (file)
index 0000000..063bc75
--- /dev/null
@@ -0,0 +1,32 @@
+AliAnalysisTaskCTauPbPb* 
+AddTaskCTau(Double_t min=0., Double_t max=90., 
+TString name="cTau_0090", Bool_t isMC=kFALSE) 
+{
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr) {
+    ::Error("AddTaskCTau", "No analysis manager to connect to.");
+    return NULL;
+  }  
+  
+  if (!mgr->GetInputEventHandler()) {
+ ::Error("AddTaskCTau","This task requires an input event handler");
+    return NULL;
+  }
+
+  AliAnalysisTaskCTauPbPb *task = new AliAnalysisTaskCTauPbPb(name);
+  task->SetCentrality(min,max);
+  task->SetMC(isMC);
+  mgr->AddTask(task);
+  
+  AliAnalysisDataContainer *cinput0 = mgr->GetCommonInputContainer();
+  mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
+
+  if (isMC) name+="_mc";
+
+  AliAnalysisDataContainer *coutput1 = 
+     mgr->CreateContainer(name, TList::Class(),
+     AliAnalysisManager::kOutputContainer, name+".root");
+  mgr->ConnectOutput(task,1,coutput1);
+
+  return task;
+}   
diff --git a/PWG2/SPECTRA/LambdaK0PbPb/AliAnalysisTaskCTauPbPb.cxx b/PWG2/SPECTRA/LambdaK0PbPb/AliAnalysisTaskCTauPbPb.cxx
new file mode 100644 (file)
index 0000000..0a4c999
--- /dev/null
@@ -0,0 +1,702 @@
+#include <TCanvas.h>
+#include <TTree.h>
+#include <TFile.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TH3F.h>
+#include <TPDGCode.h>
+#include <TDatabasePDG.h>
+#include <TParticlePDG.h>
+#include <TParticle.h>
+#include <TROOT.h>
+
+#include "AliESDEvent.h"
+#include "AliESDv0.h"
+#include "AliESDcascade.h"
+
+#include "AliCentrality.h"
+
+#include "AliMCEvent.h"
+#include "AliStack.h"
+
+#include "AliPID.h"
+#include "AliPIDResponse.h"
+
+#include "AliInputEventHandler.h"
+#include "AliAnalysisManager.h"
+
+#include "AliAnalysisTaskCTauPbPb.h"
+
+extern TROOT *gROOT;
+
+ClassImp(AliAnalysisTaskCTauPbPb)
+
+static Int_t    nbins=102;  // number of bins
+static Double_t lMin=0.0, lMax=100.;
+static Double_t pMin=0.0, pMax=10.;
+static Double_t yMax=0.5;
+
+
+//
+//  This is a little task for checking the c*tau of the strange particles 
+//
+
+AliAnalysisTaskCTauPbPb::AliAnalysisTaskCTauPbPb(const char *name) :
+AliAnalysisTaskSE(name),
+fIsMC(kFALSE),
+fCMin(0.),
+fCMax(90.),
+fOutput(0),
+fMult(0),
+fdEdx(0),
+fdEdxPid(0),
+
+fK0sM(0),
+fK0sSi(0),
+fK0sMC(0),
+fK0sAs(0),
+
+fLambdaM(0),
+fLambdaSi(0),
+fLambdaMC(0),
+fLambdaAs(0),
+
+fCPA(0),
+fDCA(0),
+
+fLambdaEff(0),
+fLambdaPt(0),
+
+fLambdaFromXi(0),
+fXiM(0),
+fXiSiP(0)
+{
+  // Constructor. Initialization of pointers
+  DefineOutput(1, TList::Class());
+}
+
+void AliAnalysisTaskCTauPbPb::UserCreateOutputObjects()
+{
+  fOutput = new TList(); 
+  fOutput->SetOwner();
+
+
+  fMult=new TH1F("fMult","Multiplicity",1100,0.,3300);
+  fMult->GetXaxis()->SetTitle("N tracks"); 
+  fOutput->Add(fMult);
+
+  fdEdx=new TH2F("fdEdx","dE/dx",50,0.2,3,50,0.,6.);
+  fOutput->Add(fdEdx);
+
+  fdEdxPid=new TH2F("fdEdxPid","dE/dx with PID",50,0.2,3,50,0.,6.);
+  fOutput->Add(fdEdxPid);
+
+  fK0sM = 
+  new TH2F("fK0sM", "Mass for K^{0}_{s}", nbins/2, 0.448, 0.548, 10,pMin,pMax);
+  fK0sM->GetXaxis()->SetTitle("Mass [GeV/c]"); 
+  fOutput->Add(fK0sM);
+
+  fK0sSi = 
+  new TH2F("fK0sSi","L_{T} vs p_{T} for K^{0}_{s}, side-band subtracted",
+  nbins,pMin,pMax,nbins,lMin,lMax);
+  fK0sSi->GetXaxis()->SetTitle("p_{T} [GeV/c]"); 
+  fK0sSi->GetYaxis()->SetTitle("L_{T} [cm]"); 
+  fOutput->Add(fK0sSi);
+
+  fK0sMC = 
+  new TH2F("fK0sMC","L_{T} vs p_{T} for K^{0}_{s}, from MC stack", 
+  nbins,pMin,pMax,nbins,lMin,lMax);
+  fK0sMC->GetXaxis()->SetTitle("p_{T} [GeV/c]"); 
+  fK0sMC->GetYaxis()->SetTitle("L_{T} [cm]"); 
+  fOutput->Add(fK0sMC);
+
+  fK0sAs = 
+  new TH2F("fK0sAs", "L_{T} vs p_{T} for K^{0}_{s}, associated", 
+  nbins,pMin,pMax,nbins,lMin,lMax);
+  fK0sAs->GetXaxis()->SetTitle("p_{T} [GeV/c]"); 
+  fK0sAs->GetYaxis()->SetTitle("L_{T} [cm]"); 
+  fOutput->Add(fK0sAs);
+
+  //----------------------
+
+  fLambdaM = 
+  new TH2F("fLambdaM","Mass for \\Lambda", nbins, 1.065, 1.165,nbins,pMin,pMax);
+  //new TH2F("fLambdaM","Mass for \\Lambda", nbins, 1.065, 1.165,10,0.1,1.1);
+  fLambdaM->GetXaxis()->SetTitle("Mass [GeV/c]"); 
+  fOutput->Add(fLambdaM);
+
+  fLambdaSi = 
+  new TH2F("fLambdaSi","L_{T} vs p_{T} for \\Lambda, side-band subtructed",
+  nbins,pMin,pMax,nbins,lMin,lMax);
+  fLambdaSi->GetXaxis()->SetTitle("p_{T} [GeV/c]"); 
+  fLambdaSi->GetYaxis()->SetTitle("L_{T} [cm]"); 
+  fOutput->Add(fLambdaSi);
+
+  fLambdaMC = 
+  new TH2F("fLambdaMC","c\\tau for \\Lambda, from MC stack", 
+  nbins,pMin,pMax,nbins,lMin,lMax);
+  fLambdaMC->GetXaxis()->SetTitle("p_{T} [GeV/c]"); 
+  fLambdaMC->GetYaxis()->SetTitle("L_{T} [cm]"); 
+  fOutput->Add(fLambdaMC);
+
+  fLambdaAs = 
+  new TH2F("fLambdaAs","c\\tau for \\Lambda, associated",
+  nbins,pMin,pMax,nbins,lMin,lMax);
+  fLambdaAs->GetXaxis()->SetTitle("p_{T} [GeV/c]"); 
+  fLambdaAs->GetYaxis()->SetTitle("L_{T} [cm]"); 
+  fOutput->Add(fLambdaAs);
+
+  fCPA=new TH1F("fCPA","Cosine of the pointing angle",30,0.9978,1.);
+  fOutput->Add(fCPA);
+  fDCA=new TH1F("fDCA","DCA between the daughters",30,0.,1.1);
+  fOutput->Add(fDCA);
+
+  fLambdaEff=fLambdaAs->ProjectionX();
+  fLambdaEff->SetName("fLambdaEff");
+  fLambdaEff->SetTitle("Efficiency for #Lambda");
+  fOutput->Add(fLambdaEff);
+
+  fLambdaPt=fLambdaAs->ProjectionX();
+  fLambdaPt->SetName("fLambdaPt");
+  fLambdaPt->SetTitle("Raw #Lambda pT spectrum");
+  fOutput->Add(fLambdaPt);
+
+  //----------------------
+
+  fLambdaFromXi=new TH3F("fLambdaFromXi","L_{T} vs p_{T} vs p_{T} of \\Xi for \\Lambda from Xi",
+  nbins,pMin,pMax,nbins,lMin,lMax,33,pMin,pMax+2);
+  fOutput->Add(fLambdaFromXi);
+
+  fXiM  = 
+  new TH2F("fXiM", "\\Xi mass distribution", 50, 1.271, 1.371,12,pMin,pMax+2);
+  fOutput->Add(fXiM);
+
+  fXiSiP  = new TH1F("fXiSiP", "Pt for \\Xi, side-band subracted",
+  33,pMin,pMax+2);
+  fOutput->Add(fXiSiP);
+
+
+  PostData(1, fOutput);
+}
+
+static Bool_t AcceptTrack(const AliESDtrack *t) {
+  if (!t->IsOn(AliESDtrack::kTPCrefit)) return kFALSE;
+  if (t->GetKinkIndex(0)>0) return kFALSE;
+
+  Float_t nCrossedRowsTPC = t->GetTPCClusterInfo(2,1); 
+  if (nCrossedRowsTPC < 70) return kFALSE;
+  Int_t findable=t->GetTPCNclsF();
+  if (findable <= 0) return kFALSE;
+  if (nCrossedRowsTPC/findable < 0.8) return kFALSE;
+
+  return kTRUE;   
+}
+
+static Bool_t AcceptV0(const AliESDv0 *v0, const AliESDEvent *esd) {
+
+  if (v0->GetOnFlyStatus()) return kFALSE;
+
+  if (v0->Pt() < pMin) return kFALSE;
+
+  Int_t nidx=TMath::Abs(v0->GetNindex());
+  AliESDtrack *ntrack=esd->GetTrack(nidx);
+  if (!AcceptTrack(ntrack)) return kFALSE;
+
+  Int_t pidx=TMath::Abs(v0->GetPindex());
+  AliESDtrack *ptrack=esd->GetTrack(pidx);
+  if (!AcceptTrack(ptrack)) return kFALSE;
+
+  Float_t xy,z0;
+  ntrack->GetImpactParameters(xy,z0);
+  if (TMath::Abs(xy)<0.1) return kFALSE;
+  ptrack->GetImpactParameters(xy,z0);
+  if (TMath::Abs(xy)<0.1) return kFALSE;
+
+  Double_t dca=v0->GetDcaV0Daughters();
+  if (dca>1.0) return kFALSE;
+  //if (dca>0.7) return kFALSE;
+  //if (dca>0.4) return kFALSE;
+
+  Double_t cpa=v0->GetV0CosineOfPointingAngle();
+  if (cpa<0.998) return kFALSE;
+  //if (cpa<0.99875) return kFALSE;
+  //if (cpa<0.9995) return kFALSE;
+
+  Double_t xx,yy,zz; v0->GetXYZ(xx,yy,zz);
+  Double_t r2=xx*xx + yy*yy;
+  if (r2<0.9*0.9) return kFALSE;
+  if (r2>100*100) return kFALSE;
+
+  return kTRUE;
+}
+
+static Bool_t AcceptCascade(const AliESDcascade *cs, const AliESDEvent *esd) {
+
+  if (cs->Pt() < pMin) return kFALSE;
+
+  Int_t bidx=TMath::Abs(cs->GetBindex());
+  AliESDtrack *btrack=esd->GetTrack(bidx);
+  if (!AcceptTrack(btrack)) return kFALSE;
+
+  Float_t xy,z0; 
+  btrack->GetImpactParameters(xy,z0);
+  if (TMath::Abs(xy)<0.03) return kFALSE;
+
+  const AliESDVertex *vtx=esd->GetPrimaryVertexSPD();
+  if (!vtx->GetStatus()) {
+     vtx=esd->GetPrimaryVertexTracks();
+     if (!vtx->GetStatus()) return kFALSE;
+  }
+  Double_t xv=vtx->GetXv(), yv=vtx->GetYv(), zv=vtx->GetZv();
+  if (cs->GetCascadeCosineOfPointingAngle(xv,yv,zv) < 0.999) return kFALSE;
+
+  if (cs->GetDcaXiDaughters() > 0.3) return kFALSE;
+
+  return kTRUE;
+}
+
+void AliAnalysisTaskCTauPbPb::UserExec(Option_t *)
+{
+
+  AliESDEvent *esd=(AliESDEvent *)InputEvent();
+
+  if (!esd) {
+    Printf("ERROR: esd not available");
+    return;
+  }
+
+  fMult->Fill(-100); //event counter  
+
+  // Physics selection
+  AliAnalysisManager *mgr= AliAnalysisManager::GetAnalysisManager();
+  AliInputEventHandler *hdr=(AliInputEventHandler*)mgr->GetInputEventHandler();
+  UInt_t maskIsSelected = hdr->IsEventSelected();
+  Bool_t isSelected = (maskIsSelected & AliVEvent::kMB);
+  if (!isSelected) return;
+
+  // Centrality selection
+  AliCentrality *cent=esd->GetCentrality();
+  if (!cent->IsEventInCentralityClass(fCMin,fCMax,"V0M")) return;
+
+  const AliESDVertex *vtx=esd->GetPrimaryVertexSPD();
+  if (!vtx->GetStatus()) {
+     vtx=esd->GetPrimaryVertexTracks();
+     if (!vtx->GetStatus()) return;
+  }
+  Double_t xv=vtx->GetXv(), yv=vtx->GetYv(), zv=vtx->GetZv();
+
+  if (TMath::Abs(zv) > 10.) return ;   
+  AliPIDResponse *pidResponse = hdr->GetPIDResponse(); 
+  //fMult->Fill(-100); //event counter  
+
+  //+++++++ MC
+  AliStack *stack = 0x0;
+  Double_t mcXv=0., mcYv=0., mcZv=0.;
+
+  if (fIsMC) {
+     AliMCEvent *mcEvent = MCEvent();
+     stack = mcEvent->Stack();
+     if (!stack) {
+        Printf("ERROR: stack not available");
+        return;
+     }
+
+     const AliVVertex *mcVtx=mcEvent->GetPrimaryVertex();
+
+     mcXv=mcVtx->GetX(); mcYv=mcVtx->GetY(); mcZv=mcVtx->GetZ();
+
+     Int_t ntrk=stack->GetNtrack(), ntrk0=ntrk;
+     while (ntrk--) {
+       TParticle *p0=stack->Particle(ntrk);
+       Int_t code=p0->GetPdgCode();
+       if (code != kK0Short)
+        if (code != kLambda0) continue;
+
+       Int_t plab=p0->GetFirstDaughter(), nlab=p0->GetLastDaughter();
+       if (nlab==plab) continue;
+       if (nlab<0) continue;
+       if (plab<0) continue;
+       if (nlab>=ntrk0) continue;
+       if (plab>=ntrk0) continue;
+       TParticle *part = stack->Particle(plab);
+       if (!part) continue;
+       TParticlePDG *partPDG = part->GetPDG();
+       if (!partPDG) continue;
+       Double_t charge=partPDG->Charge();
+       if (charge == 0.) continue;
+  
+       Double_t pt=p0->Pt();
+       if (pt<pMin) continue;
+       if (TMath::Abs(p0->Y())>yMax) continue;
+    
+       Double_t x=p0->Vx(), y=p0->Vy(), z=p0->Vz();
+       Double_t dx=mcXv-x, dy=mcYv-y, dz=mcZv-z;
+       Double_t l=TMath::Sqrt(dx*dx + dy*dy + dz*dz);
+
+       if (l > 0.01) continue; // secondary V0
+
+       x=part->Vx(); y=part->Vy();
+       dx=mcXv-x; dy=mcYv-y;
+       Double_t lt=TMath::Sqrt(dx*dx + dy*dy);
+
+       if (code == kK0Short) {
+          fK0sMC->Fill(pt,lt);
+       }
+       if (code == kLambda0) {
+          fLambdaMC->Fill(pt,lt);
+       }
+     }
+  }
+
+
+  Int_t ntrk=esd->GetNumberOfTracks();
+  Int_t mult=0;
+  Double_t nsig;
+  for (Int_t i=0; i<ntrk; i++) {
+    AliESDtrack *t=esd->GetTrack(i);
+    if (!t->IsOn(AliESDtrack::kTPCrefit)) continue;
+    Float_t xy,z0;
+    t->GetImpactParameters(xy,z0);
+    if (TMath::Abs(xy)>3.) continue;
+    if (TMath::Abs(z0)>3.) continue;
+    Double_t pt=t->Pt(),pz=t->Pz();
+    if (TMath::Abs(pz/pt)>0.8) continue;
+    mult++;
+
+    Double_t p=t->GetInnerParam()->GetP();
+    Double_t dedx=t->GetTPCsignal()/47.;
+    fdEdx->Fill(p,dedx,1);
+
+    nsig=pidResponse->NumberOfSigmasTPC(t,AliPID::kProton);
+    if (TMath::Abs(nsig) < 3.) fdEdxPid->Fill(p,dedx,1);
+
+  }
+  fMult->Fill(mult);
+
+
+  Int_t nv0 = esd->GetNumberOfV0s();
+  while (nv0--) {
+      AliESDv0 *v0=esd->GetV0(nv0);
+
+      if (!AcceptV0(v0,esd)) continue;
+
+      Int_t nidx=TMath::Abs(v0->GetNindex());
+      AliESDtrack *ntrack=esd->GetTrack(nidx);
+      Int_t pidx=TMath::Abs(v0->GetPindex());
+      AliESDtrack *ptrack=esd->GetTrack(pidx);
+
+      Double_t x,y,z; v0->GetXYZ(x,y,z);
+      Double_t dx=x-xv, dy=y-yv;
+      Double_t lt=TMath::Sqrt(dx*dx + dy*dy);
+
+      Double_t pt=v0->Pt();
+
+      Bool_t ctK=kTRUE; if (0.4977*lt/pt > 3*2.68) ctK=kFALSE;
+      Bool_t ctL=kTRUE; if (1.1157*lt/pt > 3*7.89) ctL=kFALSE;
+
+      //+++++++ MC
+      if (stack) {
+         Int_t ntrk=stack->GetNtrack();
+
+         Int_t nlab=TMath::Abs(ntrack->GetLabel());
+         Int_t plab=TMath::Abs(ptrack->GetLabel());
+
+         if (nlab<0) goto noas;      
+         if (nlab>=ntrk) goto noas;      
+         if (plab<0) goto noas;      
+         if (plab>=ntrk) goto noas;      
+
+         TParticle *np=stack->Particle(nlab);
+         TParticle *pp=stack->Particle(plab);
+         Int_t i0=pp->GetFirstMother();
+         //if (np->GetFirstMother() != i0) goto noas;
+
+         Int_t in0=np->GetFirstMother();
+         if (in0<0) goto noas;
+         if (in0>=ntrk) goto noas;
+         if (in0 != i0) { // did the negative daughter decay ?
+            TParticle *nnp=stack->Particle(in0);
+            if (nnp->GetFirstMother() != i0) goto noas;
+        }
+
+         if (i0<0) goto noas;
+         if (i0>=ntrk) goto noas;
+         TParticle *p0=stack->Particle(i0);
+
+         Int_t code=p0->GetPdgCode();
+         if (code != kK0Short)
+           if (code != kLambda0) goto noas;
+
+        if (p0->Pt()<pMin) goto noas;
+        if (TMath::Abs(p0->Y())>yMax ) goto noas;
+
+
+         Double_t dz=mcZv - p0->Vz(), dy=mcYv - p0->Vy(), dx=mcXv - p0->Vx();
+         Double_t l = TMath::Sqrt(dx*dx + dy*dy + dz*dz);
+
+         dx = mcXv - pp->Vx(); dy = mcYv - pp->Vy();
+         Double_t ltAs=TMath::Sqrt(dx*dx + dy*dy);
+         Double_t ptAs=p0->Pt();
+
+        if (l > 0.01) { // Secondary V0
+          if (code != kLambda0) goto noas;
+           Int_t nx=p0->GetFirstMother();
+           if (nx<0) goto noas;
+           if (nx>=ntrk) goto noas;
+           TParticle *xi=stack->Particle(nx);
+           Int_t xcode=xi->GetPdgCode();
+           if ( xcode != kXiMinus )
+            if( xcode != 3322 ) goto noas; 
+          fLambdaFromXi->Fill(ptAs,ltAs,xi->Pt());
+        } else {
+          if (code == kLambda0) {
+              if (ctL) fLambdaAs->Fill(ptAs,ltAs);
+           } else {
+              if (ctK)  fK0sAs->Fill(ptAs,ltAs);
+          } 
+        }
+      }
+      //++++++++
+
+  noas:
+
+      Double_t dca=v0->GetDcaV0Daughters();
+      Double_t cpa=v0->GetV0CosineOfPointingAngle();
+
+      Double_t mass=0., m=0., s=0.;
+      if (ctK)
+      if (TMath::Abs(v0->RapK0Short())<yMax) {
+         v0->ChangeMassHypothesis(kK0Short);
+
+         mass=v0->GetEffMass();
+         fK0sM->Fill(mass,pt);
+
+         m=TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass();
+         s=0.0044 + (0.008-0.0044)/(10-1)*(pt - 1.);
+         if (TMath::Abs(m-mass) < 3*s) {
+            fK0sSi->Fill(pt,lt);
+         }
+         if (TMath::Abs(m-mass + 4.5*s) < 1.5*s) {
+            fK0sSi->Fill(pt,lt,-1);
+         }
+         if (TMath::Abs(m-mass - 4.5*s) < 1.5*s) {
+            fK0sSi->Fill(pt,lt,-1);
+         }
+      }
+      
+      if (ctL)
+      if (TMath::Abs(v0->RapLambda())<yMax) {
+         Double_t p=ptrack->GetInnerParam()->GetP();
+         if (p<1.) {
+            nsig=pidResponse->NumberOfSigmasTPC(ptrack,AliPID::kProton);
+            if (TMath::Abs(nsig) > 3.) continue;
+        }
+         v0->ChangeMassHypothesis(kLambda0);
+
+         mass=v0->GetEffMass();
+         fLambdaM->Fill(mass,pt);
+
+         m=TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass();
+         //s=0.0027 + (0.004-0.0027)/(10-1)*(pt-1);
+         //s=0.0015 + (0.002-0.0015)/(2.6-1)*(pt-1);
+         s=0.0023 + (0.004-0.0023)/(6-1)*(pt-1);
+         if (TMath::Abs(m-mass) < 3*s) {
+            fLambdaSi->Fill(pt,lt);
+            fCPA->Fill(cpa,1);
+            fDCA->Fill(dca,1);
+         }
+         if (TMath::Abs(m-mass + 4.5*s) < 1.5*s) {
+            fLambdaSi->Fill(pt,lt,-1);
+            fCPA->Fill(cpa,-1);
+            fDCA->Fill(dca,-1);
+         }
+         if (TMath::Abs(m-mass - 4.5*s) < 1.5*s) {
+            fLambdaSi->Fill(pt,lt,-1);
+            fCPA->Fill(cpa,-1);
+            fDCA->Fill(dca,-1);
+         }
+      }
+  }
+
+  Double_t kine0;
+  Int_t ncs=esd->GetNumberOfCascades();
+  for (Int_t i=0; i<ncs; i++) {
+      AliESDcascade *cs=esd->GetCascade(i);
+
+      if (TMath::Abs(cs->RapXi()) > yMax) continue;
+      if (!AcceptCascade(cs,esd)) continue;
+
+      AliESDv0 *v0 = (AliESDv0*)cs;
+      if (TMath::Abs(v0->RapLambda()) > yMax) continue;
+      if (!AcceptV0(v0,esd)) continue;
+
+      Double_t pt=cs->Pt();
+
+      Int_t charge=cs->Charge();      
+      if (charge < 0) {         
+         Int_t pidx=TMath::Abs(v0->GetPindex());
+         AliESDtrack *ptrack=esd->GetTrack(pidx);
+         Double_t p=ptrack->GetInnerParam()->GetP();
+         if (p<1.) {
+            nsig=pidResponse->NumberOfSigmasTPC(ptrack,AliPID::kProton);
+            if (TMath::Abs(nsig) > 3.) continue;
+        }
+         cs->ChangeMassHypothesis(kine0,kXiMinus);
+         Double_t mass=cs->GetEffMassXi();
+        pt=cs->Pt();       
+         fXiM->Fill(mass,pt);
+         Double_t m=TDatabasePDG::Instance()->GetParticle(kXiMinus)->Mass();
+         //Double_t s=0.0037;
+         Double_t s=0.002 + (0.0032-0.002)/(6-1.5)*(pt-1.5);
+         if (TMath::Abs(m-mass) < 3*s) {
+            fXiSiP->Fill(pt);
+         }
+         if (TMath::Abs(m-mass + 4.5*s) < 1.5*s) {
+            fXiSiP->Fill(pt,-1);
+         }
+         if (TMath::Abs(m-mass - 4.5*s) < 1.5*s) {
+            fXiSiP->Fill(pt,-1);
+         }
+      }
+  }
+
+}
+
+void AliAnalysisTaskCTauPbPb::Terminate(Option_t *)
+{
+   // The Terminate() function is the last function to be called during
+   // a query. It always runs on the client, it can be used to present
+   // the results graphically or save the results to file.
+  
+  fOutput=(TList*)GetOutputData(1);
+  if (!fOutput) {
+     Printf("ERROR: fOutput not available");
+     return;
+  }
+  fMult = dynamic_cast<TH1F*>(fOutput->FindObject("fMult")) ; 
+  if (!fMult) {
+     Printf("ERROR: fMult not available");
+     return;
+  }
+
+  fdEdx = dynamic_cast<TH2F*>(fOutput->FindObject("fdEdx")) ; 
+  if (!fdEdx) {
+     Printf("ERROR: fdEdx not available");
+     return;
+  }
+
+  fdEdxPid = dynamic_cast<TH2F*>(fOutput->FindObject("fdEdxPid")) ; 
+  if (!fdEdxPid) {
+     Printf("ERROR: fdEdxPid not available");
+     return;
+  }
+
+
+  fK0sMC = dynamic_cast<TH2F*>(fOutput->FindObject("fK0sMC")) ; 
+  if (!fK0sMC) {
+     Printf("ERROR: fK0sMC not available");
+     return;
+  }
+  TH1D *k0sMcPx=fK0sMC->ProjectionX(); k0sMcPx->Sumw2();
+  fK0sAs = dynamic_cast<TH2F*>(fOutput->FindObject("fK0sAs")) ; 
+  if (!fK0sAs) {
+     Printf("ERROR: fK0sAs not available");
+     return;
+  }
+  TH1D *k0sAsPx=fK0sAs->ProjectionX(); 
+  k0sAsPx->Sumw2(); //k0sAsPx->Scale(0.69);
+
+
+
+  fLambdaFromXi = dynamic_cast<TH3F*>(fOutput->FindObject("fLambdaFromXi")) ; 
+  if (!fLambdaFromXi) {
+     Printf("ERROR: fLambdaFromXi not available");
+     return;
+  }
+  TH1D *lambdaFromXiPx=fLambdaFromXi->ProjectionX(); lambdaFromXiPx->Sumw2();
+
+
+  fLambdaMC = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaMC")) ; 
+  if (!fLambdaMC) {
+     Printf("ERROR: fLambdaMC not available");
+     return;
+  }
+  TH1D *lambdaMcPx=fLambdaMC->ProjectionX(); lambdaMcPx->Sumw2();
+
+  fLambdaAs = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaAs")) ; 
+  if (!fLambdaAs) {
+     Printf("ERROR: fLambdaAs not available");
+     return;
+  }
+  TH1D *lambdaAsPx=fLambdaAs->ProjectionX(); 
+  lambdaAsPx->Sumw2(); //lambdaAsPx->Scale(0.64);
+
+  fLambdaSi = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaSi")) ; 
+  if (!fLambdaSi) {
+     Printf("ERROR: fLambdaSi not available");
+     return;
+  }
+  TH1D *lambdaSiPx=fLambdaSi->ProjectionX(); 
+  lambdaSiPx->SetName("fLambdaPt");
+  lambdaSiPx->Sumw2();
+
+  fLambdaEff = dynamic_cast<TH1D*>(fOutput->FindObject("fLambdaEff")) ; 
+  if (!fLambdaEff) {
+     Printf("ERROR: fLambdaEff not available");
+     return;
+  }
+  fLambdaPt = dynamic_cast<TH1D*>(fOutput->FindObject("fLambdaPt")) ; 
+  if (!fLambdaPt) {
+     Printf("ERROR: fLambdaPt not available");
+     return;
+  }
+
+
+  if (!gROOT->IsBatch()) {
+
+    TCanvas *c1 = new TCanvas("c1","Mulitplicity");
+    c1->SetLogy();
+    fMult->DrawCopy() ;
+
+    new TCanvas("c2","dE/dx");
+    fdEdx->DrawCopy() ;
+
+    new TCanvas("c3","dE/dx with PID");
+    fdEdxPid->DrawCopy() ;
+
+    if (fIsMC) {
+       /*
+       TH1D effK(*k0sAsPx); effK.SetTitle("Efficiency for K0s");
+       effK.Divide(k0sAsPx,k0sMcPx,1,1,"b");
+       new TCanvas("c4","Efficiency for K0s");
+       effK.DrawCopy("E") ;
+       */
+
+       fLambdaEff->Divide(lambdaAsPx,lambdaMcPx,1,1,"b");
+       new TCanvas("c5","Efficiency for #Lambda");
+       fLambdaEff->DrawCopy("E") ;
+
+       lambdaSiPx->Add(lambdaFromXiPx,-1);
+       lambdaSiPx->Divide(fLambdaEff);
+
+       new TCanvas("c6","Corrected #Lambda pt");
+       lambdaSiPx->SetTitle("Corrected #Lambda pt");
+      *fLambdaPt = *lambdaSiPx; 
+       fLambdaPt->SetLineColor(2);
+       fLambdaPt->DrawCopy("E");
+    
+       lambdaMcPx->DrawCopy("same");
+    } else {
+       new TCanvas("c6","Raw #Lambda pt");
+       lambdaSiPx->SetTitle("Raw #Lambda pt");
+      *fLambdaPt = *lambdaSiPx; 
+       fLambdaPt->SetLineColor(2);
+       fLambdaPt->DrawCopy("E");
+    }
+  }
+}
diff --git a/PWG2/SPECTRA/LambdaK0PbPb/AliAnalysisTaskCTauPbPb.h b/PWG2/SPECTRA/LambdaK0PbPb/AliAnalysisTaskCTauPbPb.h
new file mode 100644 (file)
index 0000000..89b29a2
--- /dev/null
@@ -0,0 +1,70 @@
+#ifndef AliAnalysisTaskCTauPbPb_h
+#define AliAnalysisTaskCTauPbPb_h
+
+#include "AliAnalysisTaskSE.h"
+
+class TH1F;
+class TH2F;
+class TH3F;
+class TList;
+
+//
+//  This is a little task for checking the c*tau of the strange particles 
+//
+
+class AliAnalysisTaskCTauPbPb : public AliAnalysisTaskSE {
+
+public:
+
+  AliAnalysisTaskCTauPbPb(const char *name = "AliAnalysisTaskCTauPbPb");
+  virtual ~AliAnalysisTaskCTauPbPb() {}
+
+  void SetCentrality(Double_t min, Double_t max) {fCMin=min;fCMax=max;} 
+  void SetMC(Bool_t isMC=kTRUE) {fIsMC=isMC;} 
+  
+  virtual void   UserCreateOutputObjects();
+  virtual void   UserExec(Option_t *option);
+  virtual void   Terminate(Option_t *);  
+
+
+private: 
+
+  AliAnalysisTaskCTauPbPb(const AliAnalysisTaskCTauPbPb&);           //not implemented
+  AliAnalysisTaskCTauPbPb& operator=(const AliAnalysisTaskCTauPbPb&);//not implemented 
+
+  Bool_t fIsMC;
+  Double_t fCMin;       // Min centrality
+  Double_t fCMax;       // Max centrality
+
+  TList       *fOutput; //! The list of histograms
+
+  TH1F *fMult;       //! Track multiplicity
+  TH2F* fdEdx;       //! dEdx
+  TH2F* fdEdxPid;    //! dEdx with PID
+
+
+  TH2F* fK0sM;       //! Mass for K0s
+  TH2F* fK0sSi;      //! Side-band subtracted LvsP  for K0s 
+  TH2F* fK0sMC;      //! LvsP for the K0s from the Monte Carlo stack 
+  TH2F* fK0sAs;      //! LvsP for the K0s associated with the Monte Carlo 
+
+
+  TH2F* fLambdaM;    //! Mass for Lambdas
+  TH2F* fLambdaSi;   //! Side-band subtrated LvsP for Lambda
+  TH2F* fLambdaMC;   //! LvsP for Lambdas from the Monte Carlo stack
+  TH2F* fLambdaAs;   //! LvsP for Lambdas associated with the Monte Carlo
+
+  TH1F* fCPA;   //! cos(PA) side-band subtructed
+  TH1F* fDCA;   //! DCA daughters side-band subtructed
+
+  TH1D* fLambdaEff;  //! Efficiency for Lambda  
+  TH1D* fLambdaPt;   //! Pt spectrum for Lambda
+
+  TH3F* fLambdaFromXi;//! LvsPvsPxi for Lambdas from Xis associated with MC 
+  TH2F* fXiM;         //! Mass for Xis
+  TH1F* fXiSiP;       //! Side-band subtracted Pt for reconstructed Xi
+
+  ClassDef(AliAnalysisTaskCTauPbPb,1);
+};
+
+#endif