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
fIsMC(kFALSE),
fCMin(0.),
fCMax(90.),
+fCPA(0.9975),
+fDCA(1.0),
+fTPCcr(70.),
+fTPCcrfd(0.8),
+fDCApv(0.1),
+fRmin(0.9),
+fRmax(100.),
+
fOutput(0),
fMult(0),
+fCosPA(0),
+fDtrDCA(0),
+fTPCrows(0),
+fTPCratio(0),
+fPrimDCA(0),
+fR(0),
fdEdx(0),
fdEdxPid(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),
void AliAnalysisTaskCTauPbPbaod::UserCreateOutputObjects()
{
+ Int_t nbins=100; // number of bins
+ Double_t ltMax=100.;
+ Double_t ptMax=10.;
+
fOutput = new TList();
fOutput->SetOwner();
fMult->GetXaxis()->SetTitle("N tracks");
fOutput->Add(fMult);
+ fCosPA=new TH1F("fCosPA","Cos(PA) distribution",50,0.9975,1.0005);
+ fCosPA->GetXaxis()->SetTitle("Cos(PA)");
+ fOutput->Add(fCosPA);
+
+ fDtrDCA=new TH1F("fDtrDCA","DCA between V0 daughters",50,0.0,1.5);
+ fDtrDCA->GetXaxis()->SetTitle("DCA (rel. u.)");
+ fOutput->Add(fDtrDCA);
+
+ fTPCrows=new TH1F("fTPCrows","TPC crossed pad rows",180,0.,180.);
+ fTPCrows->GetXaxis()->SetTitle("TPC crossed pad rows");
+ fOutput->Add(fTPCrows);
+
+ fTPCratio=new TH1F("fTPCratio","TPC crossed/findable pad rows",50,0.0,1.5);
+ fTPCratio->GetXaxis()->SetTitle("TPC crossed/findable pad rows");
+ fOutput->Add(fTPCratio);
+
+ fPrimDCA=new TH1F("fPrimDCA","DCA wrt the primary vertex",50,0.0,1.5);
+ fPrimDCA->GetXaxis()->SetTitle("DCA wrt the PV (cm)");
+ fOutput->Add(fPrimDCA);
+
+ fR=new TH1F("fR","Radius of the V0 vertices",101,0.0,102);
+ fR->GetXaxis()->SetTitle("R (cm)");
+ fOutput->Add(fR);
+
fdEdx=new TH2F("fdEdx","dE/dx",50,0.2,3,50,0.,6.);
fOutput->Add(fdEdx);
fOutput->Add(fdEdxPid);
fK0sM =
- new TH2F("fK0sM", "Mass for K^{0}_{s}", nbins/2,0.448,0.548,nbins,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);
//----------------------
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","L_{T} vs p_{T} 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","L_{T} vs p_{T} 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","L_{T} vs p_{T} 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","L_{T} vs p_{T} 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);
-
- fLambdaEff=fLambdaAs->ProjectionX();
- fLambdaEff->SetName("fLambdaEff");
- fLambdaEff->SetTitle("Efficiency for #Lambda");
- fOutput->Add(fLambdaEff);
+ fLambdaBarEff=fLambdaBarAs->ProjectionX();
+ fLambdaBarEff->SetName("fLambdaBarEff");
+ fLambdaBarEff->SetTitle("Efficiency for anti-#Lambda");
+ fOutput->Add(fLambdaBarEff);
- 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,33,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);
- fOutput->Add(fLambdaFromXi);
+ 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,33,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);
PostData(1, fOutput);
}
-static Bool_t AcceptTrack(const AliAODTrack *t) {
+Bool_t AliAnalysisTaskCTauPbPbaod::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;
+ if (nCrossedRowsTPC < fTPCcr) return kFALSE;
Int_t findable=t->GetTPCNclsF();
if (findable <= 0) return kFALSE;
- if (nCrossedRowsTPC/findable < 0.8) return kFALSE;
+ if (nCrossedRowsTPC/findable < fTPCcrfd) return kFALSE;
if (TMath::Abs(t->Eta()) > 0.8) return kFALSE;
+ fTPCrows->Fill(nCrossedRowsTPC);
+ fTPCratio->Fill(nCrossedRowsTPC/findable);
+
return kTRUE;
}
-static Bool_t AcceptV0(const AliAODv0 *v1, const AliAODEvent *aod) {
-
- if (v1->GetOnFlyStatus()) return kFALSE;
-
- if (v1->Pt() < pMin) return kFALSE;
+Bool_t
+AliAnalysisTaskCTauPbPbaod::AcceptV0(const AliAODv0 *v0,const AliAODEvent *aod)
+{
+ if (v0->GetOnFlyStatus()) return kFALSE;
- const AliAODTrack *ntrack1=(AliAODTrack *)v1->GetDaughter(1);
- if (!AcceptTrack(ntrack1)) return kFALSE;
+ Double_t cpa=v0->CosPointingAngle(aod->GetPrimaryVertex());
+ if (cpa < fCPA) return kFALSE;
- const AliAODTrack *ptrack1=(AliAODTrack *)v1->GetDaughter(0);
- if (!AcceptTrack(ptrack1)) return kFALSE;
+ Double_t dca=v0->DcaV0Daughters();
+ if (dca > fDCA) 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;
+ const AliAODTrack *ntrack=(AliAODTrack *)v0->GetDaughter(1);
+ if (!AcceptTrack(ntrack)) 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;
+ const AliAODTrack *ptrack=(AliAODTrack *)v0->GetDaughter(0);
+ if (!AcceptTrack(ptrack)) 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;
+ Float_t xyn=v0->DcaNegToPrimVertex();
+ if (TMath::Abs(xyn)<fDCApv) return kFALSE;
+ Float_t xyp=v0->DcaPosToPrimVertex();
+ if (TMath::Abs(xyp)<fDCApv) return kFALSE;
- Double_t xyz[3]; v1->GetSecondaryVtx(xyz);
+ Double_t xyz[3]; v0->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;
+ if (r2<fRmin*fRmin) return kFALSE;
+ if (r2>fRmax*fRmax) return kFALSE;
+
+ fCosPA->Fill(cpa);
+ fDtrDCA->Fill(dca);
+ fPrimDCA->Fill(xyn); fPrimDCA->Fill(xyp);
+ fR->Fill(TMath::Sqrt(r2));
return kTRUE;
}
-static Bool_t AcceptCascade(const AliAODcascade *cs, const AliAODEvent *aod) {
-
- if (cs->Pt() < pMin) return kFALSE;
+Bool_t AliAnalysisTaskCTauPbPbaod::AcceptCascade(const AliAODcascade *cs, const AliAODEvent *aod) {
AliAODVertex *xi = cs->GetDecayVertexXi();
const AliAODTrack *bach=(AliAODTrack *)xi->GetDaughter(0);
static Bool_t AcceptPID(const AliPIDResponse *pidResponse,
const AliAODTrack *ptrack, const TClonesArray *stack) {
- Bool_t isProton=kTRUE;
+ const AliAODPid *pid=ptrack->GetDetPid();
+ if (!pid) return kTRUE;
+ if (pid->GetTPCmomentum() > 1.) return kTRUE;
if (stack) {
// MC PID
- Int_t ntrk=stack->GetEntriesFast();
Int_t plab=TMath::Abs(ptrack->GetLabel());
- if (plab>=0)
- if (plab<ntrk) {
- AliAODMCParticle *pp=(AliAODMCParticle*)stack->UncheckedAt(plab);
- if (pp->GetPdgCode() != kProton) isProton=kFALSE;
- }
+ AliAODMCParticle *pp=(AliAODMCParticle*)(*stack)[plab];
+ if (!pp) return kTRUE;
+ if (pp->Charge() > 0) {
+ if (pp->GetPdgCode() == kProton) return kTRUE;
+ } else {
+ if (pp->GetPdgCode() == kProtonBar) return kTRUE;
+ }
} else {
// Real PID
- const AliAODPid *pid=ptrack->GetDetPid();
- if (pid) {
- if (pid->GetTPCmomentum() < 1.) {
- Double_t nsig=pidResponse->NumberOfSigmasTPC(ptrack,AliPID::kProton);
- if (TMath::Abs(nsig) > 3.) isProton=kFALSE;
- }
- }
+ Double_t nsig=pidResponse->NumberOfSigmasTPC(ptrack,AliPID::kProton);
+ if (TMath::Abs(nsig) < 3.) return kTRUE;
}
- return isProton;
+ return kFALSE;
+}
+
+static AliAODMCParticle*
+AssociateV0(const AliAODTrack *ptrack, const AliAODTrack *ntrack,
+const TClonesArray *stack, AliAODMCParticle *&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());
+ AliAODMCParticle *n=(AliAODMCParticle*)(*stack)[nlab];
+ if (!n) return 0;
+
+ Int_t plab=TMath::Abs(ptrack->GetLabel());
+ AliAODMCParticle *p=(AliAODMCParticle*)(*stack)[plab];
+ if (!p) return 0;
+
+ Int_t imp=p->GetMother(); //V0 particle, the mother of the pos. track
+ AliAODMCParticle *p0=(AliAODMCParticle*)(*stack)[imp];
+ if (!p0) return 0;
+
+ Int_t imn=n->GetMother();
+ 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 AliAnalysisTaskCTauPbPbaod::UserExec(Option_t *)
{
+ const Double_t yMax=0.5;
+ const Double_t pMin=0.0;
+ const Double_t lMax=0.001;
AliAODEvent *aod=(AliAODEvent *)InputEvent();
mcXv=mcHdr->GetVtxX(); mcYv=mcHdr->GetVtxY(); mcZv=mcHdr->GetVtxZ();
- Int_t ntrk1=stack->GetEntriesFast(), ntrk0=ntrk1;
- while (ntrk1--) {
- AliAODMCParticle *p0=(AliAODMCParticle*)stack->UncheckedAt(ntrk1);
+ Int_t ntrk=stack->GetEntriesFast();
+ while (ntrk--) {
+ AliAODMCParticle *p0=(AliAODMCParticle*)stack->UncheckedAt(ntrk);
Int_t code=p0->GetPdgCode();
if (code != kK0Short)
if (code != kLambda0)
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);
+ AliAODMCParticle *part=(AliAODMCParticle*)(*stack)[plab];
if (!part) continue;
Double_t charge=part->Charge();
if (charge == 0.) continue;
Double_t dxmc=mcXv-x, dymc=mcYv-y, dzmc=mcZv-z;
Double_t l=TMath::Sqrt(dxmc*dxmc + dymc*dymc + dzmc*dzmc);
- if (l > 0.001) continue; // secondary V0
+ if (l > lMax) continue; // secondary V0
x=part->Xv(); y=part->Yv();
dxmc=mcXv-x; dymc=mcYv-y;
Double_t lt=TMath::Sqrt(dxmc*dxmc + dymc*dymc);
- 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;
}
}
}
- //++++++++++
+ //+++++++
- Int_t ntrk2=aod->GetNumberOfTracks();
+ Int_t ntrk1=aod->GetNumberOfTracks();
Int_t mult=0;
- for (Int_t i=0; i<ntrk2; i++) {
+ for (Int_t i=0; i<ntrk1; i++) {
AliAODTrack *t=aod->GetTrack(i);
if (t->IsMuonTrack()) continue;
if (!t->IsOn(AliAODTrack::kTPCrefit)) continue;
while (nv0--) {
AliAODv0 *v0=aod->GetV0(nv0);
+ Double_t pt=TMath::Sqrt(v0->Pt2V0());
+ if (pt < pMin) continue;
+
if (!AcceptV0(v0,aod)) continue;
const AliAODTrack *ntrack=(AliAODTrack *)v0->GetDaughter(1);
Double_t dxv=xyz[0]-xv, dyv=xyz[1]-yv;
Double_t lt=TMath::Sqrt(dxv*dxv + dyv*dyv);
- 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;
-
- Bool_t isProton =AcceptPID(pidResponse, ptrack, stack);
- Bool_t isProtonBar=AcceptPID(pidResponse, ntrack, stack);
+ if (lt/pt > 3*7.89/1.1157) continue;
- //+++++++ MC
- if (stack) {
- Int_t ntrk=stack->GetEntriesFast();
+ //--- V0 switches
+ Bool_t isK0s=kTRUE;
+ Bool_t isLambda=kTRUE;
+ Bool_t isLambdaBar=kTRUE;
- Int_t nlab=TMath::Abs(ntrack->GetLabel());
- if (nlab<0) goto noas;
- if (nlab>=ntrk) goto noas;
- AliAODMCParticle *np=(AliAODMCParticle*)stack->UncheckedAt(nlab);
+ if (0.4977*lt/pt > 3*2.68) isK0s=kFALSE;
+ if (1.1157*lt/pt > 3*7.89) isLambdaBar=isLambda=kFALSE;
- Int_t plab=TMath::Abs(ptrack->GetLabel());
- if (plab<0) goto noas;
- if (plab>=ntrk) goto noas;
- AliAODMCParticle *pp=(AliAODMCParticle*)stack->UncheckedAt(plab);
+ if (!AcceptPID(pidResponse, ptrack, stack)) isLambda=kFALSE;
+ if (!AcceptPID(pidResponse, ntrack, stack)) isLambdaBar=kFALSE;
- 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)
- if (code != kLambda0Bar) goto noas;
-
- if (p0->Pt()<pMin) goto noas;
- if (TMath::Abs(p0->Y())>yMax ) goto noas;
-
- Double_t dz0=mcZv - p0->Zv(),dy0=mcYv - p0->Yv(),dx0=mcXv - p0->Xv();
- Double_t l = TMath::Sqrt(dx0*dx0 + dy0*dy0 + dz0*dz0);
-
- dx0 = mcXv - pp->Xv(); dy0 = mcYv - pp->Yv();
- Double_t ltAs=TMath::Sqrt(dx0*dx0 + dy0*dy0);
- Double_t ptAs=p0->Pt();
-
- if (l > 0.001) { // Secondary V0
- if (code != kLambda0)
- if (code != kLambda0Bar) 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 (code==kLambda0) {
- if ( xcode != kXiMinus )
- if( xcode != 3322 ) goto noas;
- fLambdaFromXi->Fill(ptAs,ltAs,xi->Pt());
- } else if (code==kLambda0Bar) {
- if ( xcode != kXiPlusBar )
- if( xcode != -3322 ) goto noas;
- 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->DcaV0Daughters();
- Double_t cpa=v0->CosPointingAngle(aod->GetPrimaryVertex());
+ Double_t yK0s=TMath::Abs(v0->RapK0Short());
+ Double_t yLam=TMath::Abs(v0->RapLambda());
+ if (yK0s > yMax) isK0s=kFALSE;
+ if (yLam > yMax) isLambda=isLambdaBar=kFALSE;
+ //---
Double_t mass=0., m=0., s=0.;
- if (ctK)
- if (TMath::Abs(v0->RapK0Short())<yMax) {
+ if (isK0s) {
mass=v0->MassK0Short();
fK0sM->Fill(mass,pt);
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);
}
}
- if (ctL)
- if (isProton)
- if (TMath::Abs(v0->RapLambda())<yMax) {
+ if (isLambda) {
mass=v0->MassLambda();
fLambdaM->Fill(mass,pt);
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) {
mass=v0->MassAntiLambda();
fLambdaBarM->Fill(mass,pt);
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
+
+ AliAODMCParticle *mcp=0;
+ AliAODMCParticle *mc0=AssociateV0(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->Xv(), dy = mcYv - mcp->Yv();
+ Double_t ltAs=TMath::Sqrt(dx*dx + dy*dy);
+
+ Double_t dz=mcZv - mc0->Zv(); dy=mcYv - mc0->Yv(); dx=mcXv - mc0->Xv();
+ 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->GetMother();
+ AliAODMCParticle *xi=(AliAODMCParticle*)(*stack)[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;
+ }
+ }
+ //++++++
+
}
Int_t ncs=aod->GetNumberOfCascades();
for (Int_t i=0; i<ncs; i++) {
AliAODcascade *cs=aod->GetCascade(i);
+ Double_t pt=TMath::Sqrt(cs->Pt2Xi());
+ if (pt < pMin) continue;
if (TMath::Abs(cs->RapXi()) > yMax) continue;
if (!AcceptCascade(cs,aod)) continue;
if (v0->RapLambda() > yMax) continue;
if (!AcceptV0(v0,aod)) continue;
- Double_t pt=TMath::Sqrt(cs->Pt2Xi());
+ //--- Cascade switches
+ Bool_t isXiMinus=kTRUE;
+ Bool_t isXiPlusBar=kTRUE;
const AliAODTrack *ptrack=(AliAODTrack *)v0->GetDaughter(0);
- Bool_t isProton=AcceptPID(pidResponse, ptrack, stack);
+ if (!AcceptPID(pidResponse, ptrack, stack)) isXiMinus=kFALSE;
const AliAODTrack *ntrack=(AliAODTrack *)v0->GetDaughter(1);
- Bool_t isProtonBar=AcceptPID(pidResponse, ntrack, stack);
+ if (!AcceptPID(pidResponse, ntrack, stack)) isXiPlusBar=kFALSE;
Int_t charge=cs->ChargeXi();
- if (isProton)
- if (charge < 0) {
+ if (charge > 0) isXiMinus=kFALSE;
+ if (charge < 0) isXiPlusBar=kFALSE;
+ //---
+
+ if (isXiMinus) {
Double_t mass=cs->MassXi();
fXiM->Fill(mass,pt);
Double_t m=TDatabasePDG::Instance()->GetParticle(kXiMinus)->Mass();
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);
+ } else {
+ isXiMinus=kFALSE;
}
if (TMath::Abs(m-mass + 4.5*s) < 1.5*s) {
fXiSiP->Fill(pt,-1);
fXiSiP->Fill(pt,-1);
}
}
- if (isProtonBar)
- if (charge > 0) {
+
+ if (isXiPlusBar) {
Double_t mass=cs->MassXi();
fXiBarM->Fill(mass,pt);
Double_t m=TDatabasePDG::Instance()->GetParticle(kXiPlusBar)->Mass();
Double_t s=0.002 + (0.0032-0.002)/(6-1.5)*(pt-1.5);
if (TMath::Abs(m-mass) < 3*s) {
fXiBarSiP->Fill(pt);
+ } else {
+ isXiPlusBar=kFALSE;
}
if (TMath::Abs(m-mass + 4.5*s) < 1.5*s) {
fXiBarSiP->Fill(pt,-1);
fXiBarSiP->Fill(pt,-1);
}
}
+
+ if (!fIsMC) continue;
+
+ //++++++ MC
+ if (!isXiMinus)
+ if (!isXiPlusBar) continue;//check MC only for the accepted cascades
+ // Here is the future association with MC
+
}
}
+ // Lambda histograms
fLambdaFromXi = dynamic_cast<TH3F*>(fOutput->FindObject("fLambdaFromXi")) ;
if (!fLambdaFromXi) {
Printf("ERROR: fLambdaFromXi not available");
}
TH1D *lambdaFromXiPx=fLambdaFromXi->ProjectionX(); lambdaFromXiPx->Sumw2();
-
fLambdaMC = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaMC")) ;
if (!fLambdaMC) {
Printf("ERROR: fLambdaMC not available");
}
+ // 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");
effK.DrawCopy("E") ;
*/
+ //+++ Lambdas
fLambdaEff->Divide(lambdaAsPx,lambdaMcPx,1,1,"b");
new TCanvas("c5","Efficiency for #Lambda");
fLambdaEff->DrawCopy("E") ;
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");
}
}
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff