]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
A new version of UserExec()
authorbelikov <belikov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 16 Feb 2012 17:43:41 +0000 (17:43 +0000)
committerbelikov <belikov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 16 Feb 2012 17:43:41 +0000 (17:43 +0000)
PWGLF/STRANGENESS/LambdaK0PbPb/AliAnalysisTaskCTauPbPb.cxx
PWGLF/STRANGENESS/LambdaK0PbPb/AliAnalysisTaskCTauPbPb.h

index 3bafff6bac29b3b468176af09ae9aa2cd4c541d0..a58ff21f5589eed7b086ba4930682b1dffb31d86 100644 (file)
@@ -31,11 +31,6 @@ extern TROOT *gROOT;
 
 ClassImp(AliAnalysisTaskCTauPbPb)
 
-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 
@@ -61,16 +56,16 @@ fLambdaSi(0),
 fLambdaMC(0),
 fLambdaAs(0),
 
+fLambdaEff(0),
+fLambdaPt(0),
+
 fLambdaBarM(0),
 fLambdaBarSi(0),
 fLambdaBarMC(0),
 fLambdaBarAs(0),
 
-fCPA(0),
-fDCA(0),
-
-fLambdaEff(0),
-fLambdaPt(0),
+fLambdaBarEff(0),
+fLambdaBarPt(0),
 
 fLambdaFromXi(0),
 fXiM(0),
@@ -86,6 +81,10 @@ fXiBarSiP(0)
 
 void AliAnalysisTaskCTauPbPb::UserCreateOutputObjects()
 {
+  Int_t    nbins=100;  // number of bins
+  Double_t ltMax=100.;
+  Double_t ptMax=10.;
+
   fOutput = new TList(); 
   fOutput->SetOwner();
 
@@ -101,27 +100,27 @@ void AliAnalysisTaskCTauPbPb::UserCreateOutputObjects()
   fOutput->Add(fdEdxPid);
 
   fK0sM = 
-  new TH2F("fK0sM", "Mass for K^{0}_{s}", nbins/2, 0.448, 0.548, 10,pMin,pMax);
+  new TH2F("fK0sM", "Mass for K^{0}_{s}", nbins/2, 0.448, 0.548, nbins,0.,ptMax);
   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);
+  nbins,0.,ptMax,nbins,0.,ltMax);
   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);
+  nbins,0.,ptMax,nbins,0.,ltMax);
   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);
+  nbins,0.,ptMax,nbins,0.,ltMax);
   fK0sAs->GetXaxis()->SetTitle("p_{T} [GeV/c]"); 
   fK0sAs->GetYaxis()->SetTitle("L_{T} [cm]"); 
   fOutput->Add(fK0sAs);
@@ -129,100 +128,109 @@ void AliAnalysisTaskCTauPbPb::UserCreateOutputObjects()
   //----------------------
 
   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,nbins,0.,ptMax);
   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);
+  nbins,0.,ptMax,nbins,0.,ltMax);
   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);
+  nbins,0.,ptMax,nbins,0.,ltMax);
   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);
+  nbins,0.,ptMax,nbins,0.,ltMax);
   fLambdaAs->GetXaxis()->SetTitle("p_{T} [GeV/c]"); 
   fLambdaAs->GetYaxis()->SetTitle("L_{T} [cm]"); 
   fOutput->Add(fLambdaAs);
 
+  //----------------------
+
+  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);
+
+  //----------------------
 
   fLambdaBarM = 
-  new TH2F("fLambdaBarM","Mass for anti-\\Lambda", nbins, 1.065, 1.165,nbins,pMin,pMax);
+  new TH2F("fLambdaBarM","Mass for anti-\\Lambda", nbins, 1.065, 1.165,nbins,0.,ptMax);
   fLambdaBarM->GetXaxis()->SetTitle("Mass [GeV/c]"); 
   fOutput->Add(fLambdaBarM);
 
   fLambdaBarSi = 
   new TH2F("fLambdaBarSi","L_{T} vs p_{T} for anti-\\Lambda, side-band subtructed",
-  nbins,pMin,pMax,nbins,lMin,lMax);
+  nbins,0.,ptMax,nbins,0.,ltMax);
   fLambdaBarSi->GetXaxis()->SetTitle("p_{T} [GeV/c]"); 
   fLambdaBarSi->GetYaxis()->SetTitle("L_{T} [cm]"); 
   fOutput->Add(fLambdaBarSi);
 
   fLambdaBarMC = 
   new TH2F("fLambdaBarMC","c\\tau for anti-\\Lambda, from MC stack", 
-  nbins,pMin,pMax,nbins,lMin,lMax);
+  nbins,0.,ptMax,nbins,0.,ltMax);
   fLambdaBarMC->GetXaxis()->SetTitle("p_{T} [GeV/c]"); 
   fLambdaBarMC->GetYaxis()->SetTitle("L_{T} [cm]"); 
   fOutput->Add(fLambdaBarMC);
 
   fLambdaBarAs = 
   new TH2F("fLambdaBarAs","c\\tau for anti-\\Lambda, associated",
-  nbins,pMin,pMax,nbins,lMin,lMax);
+  nbins,0.,ptMax,nbins,0.,ltMax);
   fLambdaBarAs->GetXaxis()->SetTitle("p_{T} [GeV/c]"); 
   fLambdaBarAs->GetYaxis()->SetTitle("L_{T} [cm]"); 
   fOutput->Add(fLambdaBarAs);
 
 
+  //----------------------
 
-  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);
+  fLambdaBarEff=fLambdaBarAs->ProjectionX();
+  fLambdaBarEff->SetName("fLambdaBarEff");
+  fLambdaBarEff->SetTitle("Efficiency for anti-#Lambda");
+  fOutput->Add(fLambdaBarEff);
 
-  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);
+  fLambdaBarPt=fLambdaBarAs->ProjectionX();
+  fLambdaBarPt->SetName("fLambdaBarPt");
+  fLambdaBarPt->SetTitle("Raw anti-#Lambda pT spectrum");
+  fOutput->Add(fLambdaBarPt);
 
   //----------------------
 
   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);
+  nbins,0.,ptMax,nbins,0.,ltMax,33,0.,ptMax+2);
   fOutput->Add(fLambdaFromXi);
 
   fXiM  = 
-  new TH2F("fXiM", "\\Xi mass distribution", 50, 1.271, 1.371,12,pMin,pMax+2);
+  new TH2F("fXiM", "\\Xi mass distribution", 50, 1.271, 1.371,12,0.,ptMax+2);
   fOutput->Add(fXiM);
 
   fXiSiP  = new TH1F("fXiSiP", "Pt for \\Xi, side-band subracted",
-  33,pMin,pMax+2);
+  33,0.,ptMax+2);
   fOutput->Add(fXiSiP);
 
 
   fLambdaBarFromXiBar=new TH3F("fLambdaBarFromXiBar","L_{T} vs p_{T} vs p_{T} of anti-\\Xi for anti-\\Lambda from anti-\\Xi",
-  nbins,pMin,pMax,nbins,lMin,lMax,33,pMin,pMax+2);
+  nbins,0.,ptMax,nbins,0.,ltMax,33,0.,ptMax+2);
   fOutput->Add(fLambdaBarFromXiBar);
 
   fXiBarM  = 
-  new TH2F("fXiBarM", "anti-\\Xi mass distribution", 50, 1.271, 1.371,12,pMin,pMax+2);
+  new TH2F("fXiBarM", "anti-\\Xi mass distribution", 50, 1.271, 1.371,12,0.,ptMax+2);
   fOutput->Add(fXiBarM);
 
   fXiBarSiP  = new TH1F("fXiBarSiP", "Pt for anti-\\Xi, side-band subracted",
-  33,pMin,pMax+2);
+  33,0.,ptMax+2);
   fOutput->Add(fXiBarSiP);
 
 
@@ -248,8 +256,6 @@ 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;
@@ -284,8 +290,6 @@ static Bool_t AcceptV0(const AliESDv0 *v0, const AliESDEvent *esd) {
 
 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;
@@ -314,19 +318,17 @@ static Bool_t AcceptPID(const AliPIDResponse *pidResponse,
   if (!par) return kTRUE;
   if (par->GetP() > 1.) return kTRUE; 
 
+
   if (stack) {
     // MC PID
-    Int_t ntrk=stack->GetNtrack();
     Int_t plab=TMath::Abs(ptrack->GetLabel());
-    if (plab>=0)
-      if (plab<ntrk) {
-         TParticle *pp=stack->Particle(plab);
-         if (pp->GetPDG()->Charge() > 0) {
-            if (pp->GetPdgCode() == kProton)    return kTRUE;
-         } else {
-            if (pp->GetPdgCode() == kProtonBar) return kTRUE;
-        }
-      }
+    TParticle *pp=stack->Particle(plab);
+    if (!pp) return kTRUE;
+    if (pp->GetPDG()->Charge() > 0) {
+       if (pp->GetPdgCode() == kProton)    return kTRUE;
+    } else {
+       if (pp->GetPdgCode() == kProtonBar) return kTRUE;
+    }
   } else {
     // Real PID
     Double_t nsig=pidResponse->NumberOfSigmasTPC(ptrack,AliPID::kProton);
@@ -336,8 +338,51 @@ static Bool_t AcceptPID(const AliPIDResponse *pidResponse,
   return kFALSE; 
 }
 
+static TParticle*
+Associate(const AliESDtrack *ptrack,const AliESDtrack *ntrack,AliStack *stack,
+TParticle *&mcp) {
+  //
+  // Try to associate the V0 with the daughters ptrack and ntrack
+  // with the Monte Carlo
+  //
+  if (!stack) return 0;
+
+  Int_t nlab=TMath::Abs(ntrack->GetLabel());
+  TParticle *n=stack->Particle(nlab);
+  if (!n) return 0;
+
+  Int_t plab=TMath::Abs(ptrack->GetLabel());
+  TParticle *p=stack->Particle(plab);
+  if (!p) return 0;
+
+  Int_t imp=p->GetFirstMother();
+  if (imp<0) return 0;
+  if (imp>=stack->GetNtrack()) return 0;
+  TParticle *p0=stack->Particle(imp); // V0 particle, mother of the pos. track
+  if (!p0) return 0;
+
+  Int_t imn=n->GetFirstMother();
+  if (imp != imn) {  // Check decays of the daughters
+    return 0; // Fixme
+  }
+
+  Int_t code=p0->GetPdgCode();
+  if (code != kK0Short)
+     if (code != kLambda0)
+       if (code != kLambda0Bar) return 0;
+
+  mcp=p;
+
+  return p0;
+
+}
+
+
 void AliAnalysisTaskCTauPbPb::UserExec(Option_t *)
 {
+  const Double_t yMax=0.5;
+  const Double_t pMin=0.0;
+  const Double_t lMax=0.001;
 
   AliESDEvent *esd=(AliESDEvent *)InputEvent();
 
@@ -391,7 +436,7 @@ void AliAnalysisTaskCTauPbPb::UserExec(Option_t *)
 
      mcXv=mcVtx->GetX(); mcYv=mcVtx->GetY(); mcZv=mcVtx->GetZ();
 
-     Int_t ntrk=stack->GetNtrack(), ntrk0=ntrk;
+     Int_t ntrk=stack->GetNtrack();
      while (ntrk--) {
        TParticle *p0=stack->Particle(ntrk);
        Int_t code=p0->GetPdgCode();
@@ -401,10 +446,6 @@ void AliAnalysisTaskCTauPbPb::UserExec(Option_t *)
 
        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();
@@ -420,20 +461,23 @@ void AliAnalysisTaskCTauPbPb::UserExec(Option_t *)
        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.001) continue; // secondary V0
+       if (l > lMax) 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) {
+       switch (code) {
+       case kK0Short:
           fK0sMC->Fill(pt,lt);
-       }
-       if (code == kLambda0) {
+          break;
+       case kLambda0:
           fLambdaMC->Fill(pt,lt);
-       }
-       if (code == kLambda0Bar) {
+          break;
+       case kLambda0Bar:
           fLambdaBarMC->Fill(pt,lt);
+          break;
+       default: break;
        }
      }
   }
@@ -466,8 +510,16 @@ void AliAnalysisTaskCTauPbPb::UserExec(Option_t *)
 
   Int_t nv0 = esd->GetNumberOfV0s();
   while (nv0--) {
+
+      Bool_t isK0s=kTRUE;
+      Bool_t isLambda=kTRUE;
+      Bool_t isLambdaBar=kTRUE;
+
       AliESDv0 *v0=esd->GetV0(nv0);
 
+      Double_t pt=v0->Pt();
+      if (pt < pMin) continue;
+
       if (!AcceptV0(v0,esd)) continue;
 
       Int_t nidx=TMath::Abs(v0->GetNindex());
@@ -479,97 +531,22 @@ void AliAnalysisTaskCTauPbPb::UserExec(Option_t *)
       Double_t dx1=x-xv, dy1=y-yv;
       Double_t lt=TMath::Sqrt(dx1*dx1 + dy1*dy1);
 
-      Double_t pt=v0->Pt();
+      if (lt/pt > 3*7.89/1.1157) continue;  
 
-      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;
+      if (0.4977*lt/pt > 3*2.68) isK0s=kFALSE;
+      if (1.1157*lt/pt > 3*7.89) isLambdaBar=isLambda=kFALSE;
 
-      Bool_t isProton   =AcceptPID(pidResponse, ptrack, stack);
-      Bool_t isProtonBar=AcceptPID(pidResponse, ntrack, stack);
+      isLambda    = isLambda && AcceptPID(pidResponse, ptrack, stack);
+      isLambdaBar = isLambdaBar && AcceptPID(pidResponse, ntrack, stack);
 
-      //+++++++ MC
-      if (stack) {
-         Int_t ntrk=stack->GetNtrack();
-
-         Int_t nlab=TMath::Abs(ntrack->GetLabel());
-         if (nlab<0) goto noas;      
-         if (nlab>=ntrk) goto noas;      
-         TParticle *np=stack->Particle(nlab);
-
-         Int_t plab=TMath::Abs(ptrack->GetLabel());
-         if (plab<0) goto noas;      
-         if (plab>=ntrk) goto noas;      
-         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)
-              if (code != kLambda0Bar) 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.001) { // Secondary V0
-          if (code != kLambda0)
-            if (code != kLambda0Bar) 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 (code == kLambda0) {
-              if ( xcode != kXiMinus )
-                if ( xcode != 3322 ) goto noas; 
-             if (ctL) fLambdaFromXi->Fill(ptAs,ltAs,xi->Pt());
-          } else if (code == kLambda0Bar) {
-              if ( xcode != kXiPlusBar )
-                if ( xcode != -3322 ) goto noas; 
-             if (ctL) fLambdaBarFromXiBar->Fill(ptAs,ltAs,xi->Pt());
-          }
-        } else {
-          if (code == kLambda0) {
-              if (ctL) fLambdaAs->Fill(ptAs,ltAs);
-          } else if (code == kLambda0Bar) {
-              if (ctL) fLambdaBarAs->Fill(ptAs,ltAs);
-           } else {
-              if (ctK)  fK0sAs->Fill(ptAs,ltAs);
-          } 
-        }
-      }
-      //++++++++
-
-  noas:
-
-      Double_t dca=v0->GetDcaV0Daughters();
-      Double_t cpa=v0->GetV0CosineOfPointingAngle();
+      Double_t yK0s=TMath::Abs(v0->RapK0Short());
+      Double_t yLam=TMath::Abs(v0->RapLambda());
+      isK0s       = isK0s && (yK0s < yMax);
+      isLambda    = isLambda && (yLam < yMax);
+      isLambdaBar = isLambdaBar && (yLam < yMax);
 
       Double_t mass=0., m=0., s=0.;
-      if (ctK)
-      if (TMath::Abs(v0->RapK0Short())<yMax) {
+      if (isK0s) {
          v0->ChangeMassHypothesis(kK0Short);
 
          mass=v0->GetEffMass();
@@ -579,6 +556,8 @@ void AliAnalysisTaskCTauPbPb::UserExec(Option_t *)
          s=0.0044 + (0.008-0.0044)/(10-1)*(pt - 1.);
          if (TMath::Abs(m-mass) < 3*s) {
             fK0sSi->Fill(pt,lt);
+         } else {
+           isK0s=kFALSE;
          }
          if (TMath::Abs(m-mass + 4.5*s) < 1.5*s) {
             fK0sSi->Fill(pt,lt,-1);
@@ -588,9 +567,7 @@ void AliAnalysisTaskCTauPbPb::UserExec(Option_t *)
          }
       }
       
-      if (ctL)
-      if (isProton)
-      if (TMath::Abs(v0->RapLambda())<yMax) {
+      if (isLambda) {
          v0->ChangeMassHypothesis(kLambda0);
 
          mass=v0->GetEffMass();
@@ -602,24 +579,18 @@ void AliAnalysisTaskCTauPbPb::UserExec(Option_t *)
          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);
-         }
+         } else {
+           isLambda=kFALSE;
+         } 
          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);
          }
       }
 
-      if (ctL)
-      if (isProtonBar)
-      if (TMath::Abs(v0->RapLambda())<yMax) {
+      if (isLambdaBar) {
          v0->ChangeMassHypothesis(kLambda0Bar);
 
          mass=v0->GetEffMass();
@@ -631,20 +602,79 @@ void AliAnalysisTaskCTauPbPb::UserExec(Option_t *)
          s=0.0023 + (0.004-0.0023)/(6-1)*(pt-1);
          if (TMath::Abs(m-mass) < 3*s) {
             fLambdaBarSi->Fill(pt,lt);
-            fCPA->Fill(cpa,1);
-            fDCA->Fill(dca,1);
+         } else {
+           isLambdaBar=kFALSE;
          }
          if (TMath::Abs(m-mass + 4.5*s) < 1.5*s) {
             fLambdaBarSi->Fill(pt,lt,-1);
-            fCPA->Fill(cpa,-1);
-            fDCA->Fill(dca,-1);
          }
          if (TMath::Abs(m-mass - 4.5*s) < 1.5*s) {
             fLambdaBarSi->Fill(pt,lt,-1);
-            fCPA->Fill(cpa,-1);
-            fDCA->Fill(dca,-1);
          }
       }
+
+      if (!fIsMC) continue;
+
+      //++++++ MC 
+      if (!isK0s)
+         if (!isLambda)
+             if (!isLambdaBar) continue;//check MC only for the accepted V0s 
+
+      TParticle *mcp=0;
+      TParticle *mc0=Associate(ptrack,ntrack,stack,mcp);
+      if (!mc0) continue;
+
+      Double_t ptAs=mc0->Pt();
+      if (ptAs < pMin) continue;
+      Double_t yAs=mc0->Y();
+      if (TMath::Abs(yAs) > yMax) continue;
+
+      Int_t code=mc0->GetPdgCode();
+
+      Double_t dx = mcXv - mcp->Vx(), dy = mcYv - mcp->Vy();
+      Double_t ltAs=TMath::Sqrt(dx*dx + dy*dy);
+      Double_t dz=mcZv - mc0->Vz(); dy=mcYv - mc0->Vy(); dx=mcXv - mc0->Vx();
+      Double_t l = TMath::Sqrt(dx*dx + dy*dy + dz*dz);
+      if (l<lMax) { // Primary V0s
+        switch (code) {
+         case kK0Short:
+            if (isK0s)       fK0sAs->Fill(ptAs,ltAs);
+            break;
+        case kLambda0:
+            if (isLambda)    fLambdaAs->Fill(ptAs,ltAs);
+            break;
+        case kLambda0Bar:
+            if (isLambdaBar) fLambdaBarAs->Fill(ptAs,ltAs);
+            break;
+         default: break;
+        }
+      } else {
+        if (code==kK0Short) continue;
+
+         Int_t nx=mc0->GetFirstMother();
+         if (nx<0) continue;
+         if (nx>=stack->GetNtrack()) continue;
+         TParticle *xi=stack->Particle(nx);
+         if (!xi) continue;
+         Int_t xcode=xi->GetPdgCode();
+         
+        switch (code) {
+        case kLambda0:
+            if (isLambda)
+           if ((xcode==kXiMinus) || (xcode==3322)) 
+               fLambdaFromXi->Fill(ptAs,ltAs,xi->Pt());
+            break;
+        case kLambda0Bar:
+            if (isLambdaBar)
+           if ((xcode==kXiPlusBar)||(xcode==-3322)) 
+               fLambdaBarFromXiBar->Fill(ptAs,ltAs,xi->Pt());
+            break;
+         default: break;
+        }
+      }
+      //++++++
+  
   }
 
   Double_t kine0;
@@ -652,10 +682,12 @@ void AliAnalysisTaskCTauPbPb::UserExec(Option_t *)
   for (Int_t i=0; i<ncs; i++) {
       AliESDcascade *cs=esd->GetCascade(i);
 
+      if (cs->Pt() < pMin) continue;
       if (TMath::Abs(cs->RapXi()) > yMax) continue;
       if (!AcceptCascade(cs,esd)) continue;
 
       AliESDv0 *v0 = (AliESDv0*)cs;
+      if (v0->Pt() < pMin) continue;
       if (TMath::Abs(v0->RapLambda()) > yMax) continue;
       if (!AcceptV0(v0,esd)) continue;
 
@@ -759,6 +791,7 @@ void AliAnalysisTaskCTauPbPb::Terminate(Option_t *)
 
 
 
+  // Lambda histograms 
   fLambdaFromXi = dynamic_cast<TH3F*>(fOutput->FindObject("fLambdaFromXi")) ; 
   if (!fLambdaFromXi) {
      Printf("ERROR: fLambdaFromXi not available");
@@ -766,7 +799,6 @@ void AliAnalysisTaskCTauPbPb::Terminate(Option_t *)
   }
   TH1D *lambdaFromXiPx=fLambdaFromXi->ProjectionX(); lambdaFromXiPx->Sumw2();
 
-
   fLambdaMC = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaMC")) ; 
   if (!fLambdaMC) {
      Printf("ERROR: fLambdaMC not available");
@@ -803,6 +835,52 @@ void AliAnalysisTaskCTauPbPb::Terminate(Option_t *)
   }
 
 
+  // anti-Lambda histograms 
+  fLambdaBarFromXiBar = 
+  dynamic_cast<TH3F*>(fOutput->FindObject("fLambdaBarFromXiBar")) ; 
+  if (!fLambdaBarFromXiBar) {
+     Printf("ERROR: fLambdaBarFromXiBar not available");
+     return;
+  }
+  TH1D *lambdaBarFromXiBarPx=
+  fLambdaBarFromXiBar->ProjectionX("Bar"); lambdaBarFromXiBarPx->Sumw2();
+
+  fLambdaBarMC = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaBarMC")) ; 
+  if (!fLambdaBarMC) {
+     Printf("ERROR: fLambdaBarMC not available");
+     return;
+  }
+  TH1D *lambdaBarMcPx=fLambdaBarMC->ProjectionX(); lambdaBarMcPx->Sumw2();
+
+  fLambdaBarAs = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaBarAs")) ; 
+  if (!fLambdaBarAs) {
+     Printf("ERROR: fLambdaBarAs not available");
+     return;
+  }
+  TH1D *lambdaBarAsPx=fLambdaBarAs->ProjectionX(); 
+  lambdaBarAsPx->Sumw2(); //lambdaBarAsPx->Scale(0.64);
+
+  fLambdaBarSi = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaBarSi")) ; 
+  if (!fLambdaBarSi) {
+     Printf("ERROR: fLambdaBarSi not available");
+     return;
+  }
+  TH1D *lambdaBarSiPx=fLambdaBarSi->ProjectionX(); 
+  lambdaBarSiPx->SetName("fLambdaBarPt");
+  lambdaBarSiPx->Sumw2();
+
+  fLambdaBarEff = dynamic_cast<TH1D*>(fOutput->FindObject("fLambdaBarEff")) ; 
+  if (!fLambdaBarEff) {
+     Printf("ERROR: fLambdaBarEff not available");
+     return;
+  }
+  fLambdaBarPt = dynamic_cast<TH1D*>(fOutput->FindObject("fLambdaBarPt")) ; 
+  if (!fLambdaBarPt) {
+     Printf("ERROR: fLambdaBarPt not available");
+     return;
+  }
+
+
   if (!gROOT->IsBatch()) {
 
     TCanvas *c1 = new TCanvas("c1","Mulitplicity");
@@ -823,6 +901,7 @@ void AliAnalysisTaskCTauPbPb::Terminate(Option_t *)
        effK.DrawCopy("E") ;
        */
 
+       //+++ Lambdas
        fLambdaEff->Divide(lambdaAsPx,lambdaMcPx,1,1,"b");
        new TCanvas("c5","Efficiency for #Lambda");
        fLambdaEff->DrawCopy("E") ;
@@ -838,12 +917,35 @@ void AliAnalysisTaskCTauPbPb::Terminate(Option_t *)
     
        lambdaMcPx->DrawCopy("same");
  
+
+       //+++ anti-Lambdas
+       fLambdaBarEff->Divide(lambdaBarAsPx,lambdaBarMcPx,1,1,"b");
+       new TCanvas("c7","Efficiency for anti-#Lambda");
+       fLambdaBarEff->DrawCopy("E") ;
+
+       lambdaBarSiPx->Add(lambdaBarFromXiBarPx,-1);
+       lambdaBarSiPx->Divide(fLambdaBarEff);
+
+       new TCanvas("c8","Corrected anti-#Lambda pt");
+       lambdaBarSiPx->SetTitle("Corrected anti-#Lambda pt");
+      *fLambdaBarPt = *lambdaBarSiPx; 
+       fLambdaBarPt->SetLineColor(2);
+       fLambdaBarPt->DrawCopy("E");
+    
+       lambdaBarMcPx->DrawCopy("same");
     } else {
        new TCanvas("c6","Raw #Lambda pt");
        lambdaSiPx->SetTitle("Raw #Lambda pt");
       *fLambdaPt = *lambdaSiPx; 
        fLambdaPt->SetLineColor(2);
        fLambdaPt->DrawCopy("E");
+
+
+       new TCanvas("c7","Raw anti-#Lambda pt");
+       lambdaBarSiPx->SetTitle("Raw anti-#Lambda pt");
+      *fLambdaBarPt = *lambdaBarSiPx; 
+       fLambdaBarPt->SetLineColor(2);
+       fLambdaBarPt->DrawCopy("E");
     }
   }
 }
index 6f7a67653841fc05da0d905f995d41a676776a6d..454119da7a35d4015cccc0c080dbd5091c05d5db 100644 (file)
@@ -54,16 +54,16 @@ private:
   TH2F* fLambdaMC;   //! LvsP for Lambdas from the Monte Carlo stack
   TH2F* fLambdaAs;   //! LvsP for Lambdas associated with the Monte Carlo
 
+  TH1D* fLambdaEff;  //! Efficiency for Lambda  
+  TH1D* fLambdaPt;   //! Pt spectrum for Lambda
+
   TH2F* fLambdaBarM;  //! Mass for anti-Lambdas
   TH2F* fLambdaBarSi; //! Side-band subtrated LvsP for anti-Lambda
   TH2F* fLambdaBarMC; //! LvsP for anti-Lambdas from the Monte Carlo stack
   TH2F* fLambdaBarAs; //! LvsP for anti-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
+  TH1D* fLambdaBarEff;  //! Efficiency for anti-Lambda  
+  TH1D* fLambdaBarPt;   //! Pt spectrum for anti-Lambda
 
   TH3F* fLambdaFromXi;//! LvsPvsPxi for Lambdas from Xis associated with MC 
   TH2F* fXiM;         //! Mass for Xis
@@ -73,7 +73,7 @@ private:
   TH2F* fXiBarM;         //! Mass for anti-Xis
   TH1F* fXiBarSiP;       //! Side-band subtracted Pt for reconstructed anti-Xi
 
-  ClassDef(AliAnalysisTaskCTauPbPb,2);
+  ClassDef(AliAnalysisTaskCTauPbPb,3);
 };
 
 #endif