[u/mrichter/AliRoot.git] / PWG2 / SPECTRA / LambdaK0PbPb / AliAnalysisTaskCTauPbPbaod.cxx

#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 <TClonesArray.h>
#include <TROOT.h>

#include "AliAODMCHeader.h"
#include "AliAODMCParticle.h"

#include "AliAODEvent.h"
#include "AliAODv0.h"
#include "AliAODcascade.h"

#include "AliCentrality.h"

#include "AliPID.h"
#include "AliPIDResponse.h"
#include "AliAODPid.h"

#include "AliInputEventHandler.h"
#include "AliAnalysisManager.h"

#include "AliAnalysisTaskCTauPbPbaod.h"

extern TROOT *gROOT;

ClassImp(AliAnalysisTaskCTauPbPbaod)

static Int_t    nbins=100;  // 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 
//

AliAnalysisTaskCTauPbPbaod::AliAnalysisTaskCTauPbPbaod(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 AliAnalysisTaskCTauPbPbaod::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,nbins,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);
  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","L_{T} vs p_{T} 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","L_{T} vs p_{T} 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 AliAODTrack *t) {
  if (!t->IsOn(AliAODTrack::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 AliAODv0 *v1, const AliAODEvent *aod) {

  if (v1->GetOnFlyStatus()) return kFALSE;

  if (v1->Pt() < pMin) return kFALSE;

  const AliAODTrack *ntrack1=(AliAODTrack *)v1->GetDaughter(1);
  if (!AcceptTrack(ntrack1)) return kFALSE;

  const AliAODTrack *ptrack1=(AliAODTrack *)v1->GetDaughter(0);
  if (!AcceptTrack(ptrack1)) return kFALSE;

  Float_t xy=v1->DcaNegToPrimVertex();
  if (TMath::Abs(xy)<0.1) return kFALSE;
  xy=v1->DcaPosToPrimVertex();
  if (TMath::Abs(xy)<0.1) return kFALSE;

  Double_t dca=v1->DcaV0Daughters();
  if (dca>1.0) return kFALSE;
  //if (dca>0.7) return kFALSE;
  //if (dca>0.4) return kFALSE;

  Double_t cpa=v1->CosPointingAngle(aod->GetPrimaryVertex());
  if (cpa<0.998) return kFALSE;
  //if (cpa<0.99875) return kFALSE;
  //if (cpa<0.9995) return kFALSE;

  Double_t xyz[3]; v1->GetSecondaryVtx(xyz);
  Double_t r2=xyz[0]*xyz[0] + xyz[1]*xyz[1];
  if (r2<0.9*0.9) return kFALSE;
  if (r2>100*100) return kFALSE;

  return kTRUE;
}

static Bool_t AcceptCascade(const AliAODcascade *cs, const AliAODEvent *aod) {

  if (cs->Pt() < pMin) return kFALSE;

  AliAODVertex *xi = cs->GetDecayVertexXi(); 
  const AliAODTrack *bach=(AliAODTrack *)xi->GetDaughter(0);
  if (!AcceptTrack(bach)) return kFALSE;
   
  Double_t xy=cs->DcaBachToPrimVertex();
  if (TMath::Abs(xy) < 0.03) return kFALSE;

  const AliAODVertex *vtx=aod->GetPrimaryVertex();
  Double_t xv=vtx->GetX(), yv=vtx->GetY(), zv=vtx->GetZ();
  Double_t cpa=cs->CosPointingAngleXi(xv,yv,zv);
  if (cpa<0.999) return kFALSE;

  if (cs->DcaXiDaughters() > 0.3) return kFALSE;

  return kTRUE;
}

void AliAnalysisTaskCTauPbPbaod::UserExec(Option_t *)
{

  AliAODEvent *aod=(AliAODEvent *)InputEvent();

  if (!aod) {
    Printf("ERROR: aod 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=aod->GetCentrality();
  if (!cent->IsEventInCentralityClass(fCMin,fCMax,"V0M")) return;

  const AliAODVertex *vtx=aod->GetPrimaryVertex();
  if (vtx->GetNContributors()<3) return;

  Double_t xv=vtx->GetX(), yv=vtx->GetY(), zv=vtx->GetZ();

  if (TMath::Abs(zv) > 10.) return ;   
 
  AliPIDResponse *pidResponse = hdr->GetPIDResponse(); 
 
  //fMult->Fill(-100); //event counter  

  //+++++++ MC
  TClonesArray *stack = 0x0;
  Double_t mcXv=0., mcYv=0., mcZv=0.;
   
  if (fIsMC) {
     TList *lst = aod->GetList();
     stack = (TClonesArray*)lst->FindObject(AliAODMCParticle::StdBranchName());
     if (!stack) {
        Printf("ERROR: stack not available");
        return;
     }
     AliAODMCHeader *
     mcHdr=(AliAODMCHeader*)lst->FindObject(AliAODMCHeader::StdBranchName());

     mcXv=mcHdr->GetVtxX(); mcYv=mcHdr->GetVtxY(); mcZv=mcHdr->GetVtxZ();

     Int_t ntrk=stack->GetEntriesFast(), ntrk0=ntrk;
     while (ntrk--) {
       AliAODMCParticle *p0=(AliAODMCParticle*)stack->UncheckedAt(ntrk);
       Int_t code=p0->GetPdgCode();
       if (code != kK0Short)
	 if (code != kLambda0) continue;

       Int_t plab=p0->GetDaughter(0), nlab=p0->GetDaughter(1);
       if (nlab==plab) continue;
       if (nlab<0) continue;
       if (plab<0) continue;
       if (nlab>=ntrk0) continue;
       if (plab>=ntrk0) continue;
       AliAODMCParticle *part = (AliAODMCParticle *)stack->UncheckedAt(plab);
       if (!part) continue;
       Double_t charge=part->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->Xv(), y=p0->Yv(), z=p0->Zv();
       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->Xv(); y=part->Yv();
       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=aod->GetNumberOfTracks();
  Int_t mult=0;
  Double_t nsig;
  for (Int_t i=0; i<ntrk; i++) {
    AliAODTrack *t=aod->GetTrack(i);
    if (t->IsMuonTrack()) continue;
    if (!t->IsOn(AliAODTrack::kTPCrefit)) continue;

    Double_t xyz[3];
    if (t->GetPosition(xyz)) continue;
    if (TMath::Abs(xyz[0])>3.) continue;
    if (TMath::Abs(xyz[1])>3.) continue;

    Double_t pt=t->Pt(),pz=t->Pz();
    if (TMath::Abs(pz/pt)>0.8) continue;

    mult++;

    const AliAODPid *pid=t->GetDetPid();
    if (!pid) continue;

    Double_t p=pid->GetTPCmomentum();
    Double_t dedx=pid->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 = aod->GetNumberOfV0s();
  while (nv0--) {
      AliAODv0 *v0=aod->GetV0(nv0);

      if (!AcceptV0(v0,aod)) continue;
      
      const AliAODTrack *ntrack=(AliAODTrack *)v0->GetDaughter(1);
      const AliAODTrack *ptrack=(AliAODTrack *)v0->GetDaughter(0);

      Double_t xyz[3]; v0->GetSecondaryVtx(xyz);
      Double_t dx=xyz[0]-xv, dy=xyz[1]-yv;
      Double_t lt=TMath::Sqrt(dx*dx + dy*dy);

      Double_t pt=TMath::Sqrt(v0->Pt2V0());

      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->GetEntriesFast();

         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;      

         AliAODMCParticle *np=(AliAODMCParticle*)stack->UncheckedAt(nlab);
         AliAODMCParticle *pp=(AliAODMCParticle*)stack->UncheckedAt(plab);

         Int_t i0=pp->GetMother();
         //if (np->GetMother() != i0) goto noas;

         Int_t in0=np->GetMother();
         if (in0<0) goto noas;
         if (in0>=ntrk) goto noas;
         if (in0 != i0) { // did the negative daughter decay ?
            AliAODMCParticle *nnp=(AliAODMCParticle*)stack->UncheckedAt(in0);
            if (nnp->GetMother() != i0) goto noas;
	 }

         if (i0<0) goto noas;
         if (i0>=ntrk) goto noas;
         AliAODMCParticle *p0=(AliAODMCParticle*)stack->UncheckedAt(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->Zv(), dy=mcYv - p0->Yv(), dx=mcXv - p0->Xv();
         Double_t l = TMath::Sqrt(dx*dx + dy*dy + dz*dz);

         dx = mcXv - pp->Xv(); dy = mcYv - pp->Yv();
         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->GetMother();
           if (nx<0) goto noas;
           if (nx>=ntrk) goto noas;
           AliAODMCParticle *xi=(AliAODMCParticle*)stack->UncheckedAt(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->DcaV0Daughters();
      Double_t cpa=v0->CosPointingAngle(aod->GetPrimaryVertex());

      Double_t mass=0., m=0., s=0.;
      if (ctK)
      if (TMath::Abs(v0->RapK0Short())<yMax) {
         mass=v0->MassK0Short();
         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) {
         const AliAODPid *pid=ptrack->GetDetPid();
         if (pid) {
            Double_t p=pid->GetTPCmomentum();
            if (p<1.) {
               nsig=pidResponse->NumberOfSigmasTPC(ptrack,AliPID::kProton);
               if (TMath::Abs(nsig) > 3.) continue;
	    }
	 }
         mass=v0->MassLambda();
         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);
         }
      }
  }

  Int_t ncs=aod->GetNumberOfCascades();
  for (Int_t i=0; i<ncs; i++) {
      AliAODcascade *cs=aod->GetCascade(i);
 
      if (TMath::Abs(cs->RapXi()) > yMax) continue;
      if (!AcceptCascade(cs,aod)) continue;

      const AliAODv0 *v0=(AliAODv0*)cs;
      if (v0->RapLambda() > yMax) continue;
      if (!AcceptV0(v0,aod)) continue;

      Double_t pt=TMath::Sqrt(cs->Pt2Xi());

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