//********************************************************************
// Example (very naive for the moment) of the data analysis
// using the ESD classes
-// It demonstrates the idea of the "combined PID".
-// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
+// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
//********************************************************************
-#ifndef __CINT__
+#if !defined( __CINT__) || defined(__MAKECINT__)
#include <Riostream.h>
- #include "TKey.h"
+ #include <TTree.h>
#include "TFile.h"
- #include "TH2F.h"
+ #include "TH1F.h"
#include "TCanvas.h"
+ #include "TStyle.h"
#include "TStopwatch.h"
- #include "AliESD.h"
+ #include "AliESDEvent.h"
#endif
-Int_t AliESDanalysis(Int_t nev=1) {
- TH2F *tpcHist=
- new TH2F("tpcHist","TPC dE/dX vs momentum",100,0.,4.,100,0.,500.);
- tpcHist->SetXTitle("p (GeV/c)"); tpcHist->SetYTitle("dE/dx (Arb. Units)");
- tpcHist->SetMarkerStyle(8);
- tpcHist->SetMarkerSize(0.3);
-
- TH2F *elHist=new TH2F("elHist","dE/dX vs momentum",100,0.,4.,50,0.,500.);
- elHist->SetMarkerStyle(8);
- elHist->SetMarkerSize(0.3);
-
- TH2F *piHist=new TH2F("piHist","dE/dX vs momentum",100,0.,4.,50,0.,500.);
- piHist->SetMarkerColor(2);
- piHist->SetMarkerStyle(8);
- piHist->SetMarkerSize(0.3);
-
- TH2F *kaHist=new TH2F("kaHist","dE/dX vs momentum",100,0.,4.,100,0.,500.);
- kaHist->SetMarkerColor(4);
- kaHist->SetMarkerStyle(8);
- kaHist->SetMarkerSize(0.3);
-
- TH2F *prHist=new
- TH2F("prHist","Classification into e, pi, K and p",100,0.,4.,100,0.,500.);
- prHist->SetXTitle("p (GeV/c)"); prHist->SetYTitle("dE/dx (Arb. Units)");
- prHist->SetMarkerColor(6);
- prHist->SetMarkerStyle(8);
- prHist->SetMarkerSize(0.3);
-
- TFile *ef=TFile::Open("AliESDs.root");
- if (!ef->IsOpen()) {cerr<<"Can't AliESDs.root !\n"; return 1;}
+extern TStyle *gStyle;
+Int_t AliESDanalysis() {
TStopwatch timer;
- Int_t rc=0,n=0;
- TKey *key=0;
- TIter next(ef->GetListOfKeys());
- //****** Tentative particle type "concentrations"
- Double_t c[5]={0.05, 0., 0.85, 0.10, 0.05};
+ gStyle->SetOptStat(111110);
+ gStyle->SetOptFit(1);
+
+ Double_t V0mass=0.497672, V0width=0.020, V0window=0.05;
+ Double_t mmin=V0mass-V0window, mmax=V0mass+V0window;
+ TH1F *hm =new TH1F("hm","K0s",40, mmin, mmax);
+ hm->SetXTitle("Mass (GeV/c**2)"); hm->SetLineColor(2);
+ TH1F *hp =new TH1F("hp","Momentum of the positive daughter",40, 0, 2);
+ hp->SetXTitle("P (GeV/c)"); hp->SetLineColor(4);
+
+//****** File with the ESD
+ TFile *ef=TFile::Open("AliESDs.root");
+ if (!ef || !ef->IsOpen()) {cerr<<"Can't AliESDs.root !\n"; return 1;}
+ AliESDEvent* event = new AliESDEvent();
+ TTree* tree = (TTree*) ef->Get("esdTree");
+ if (!tree) {cerr<<"no ESD tree found\n"; return 1;};
+ event->ReadFromTree(tree);
- //******* The loop over events
- while ((key=(TKey*)next())!=0) {
- cerr<<"Processing event number : "<<n++<<endl;
- AliESD *event=(AliESD*)key->ReadObj();
+ Int_t n=0;
+
+//******* The loop over events
+ while (tree->GetEvent(n)) {
+ cout<<endl<<"Processing event number : "<<n++<<endl;
Int_t ntrk=event->GetNumberOfTracks();
- cerr<<"Number of ESD tracks : "<<ntrk<<endl;
+ cout<<"Number of ESD tracks : "<<ntrk<<endl;
+ Int_t nv0=event->GetNumberOfV0s();
+ cout<<"Number of ESD V0s : "<<nv0<<endl;
+ Int_t ncas=event->GetNumberOfCascades();
+ cout<<"Number of ESD cascades : "<<ncas<<endl;
+
//****** The loop over tracks
+ Int_t nk=0;
while (ntrk--) {
- AliESDtrack *t=event->GetTrack(ntrk);
-
- Double_t p=t->GetP();
-
- if (t->GetStatus()&AliESDtrack::kTPCin) {
- Double_t dedx=t->GetTPCsignal();
- tpcHist->Fill(p,dedx,1);
- }
-
- if (t->GetStatus()&AliESDtrack::kESDpid) {
- Double_t dedx=t->GetTPCsignal();
- Double_t r[10]; t->GetESDpid(r);
- Double_t rc=0.;
- Int_t i;
- for (i=0; i<AliESDtrack::kSPECIES; i++) rc+=(c[i]*r[i]);
- if (rc==0.) continue;
-
- //Here we apply Bayes' formula
- Double_t w[10];
- for (i=0; i<AliESDtrack::kSPECIES; i++) w[i]=c[i]*r[i]/rc;
-
- if (w[4]>w[3] && w[4]>w[2] && w[4]>w[0]) prHist->Fill(p,dedx,1);
- if (w[3]>w[4] && w[3]>w[2] && w[3]>w[0]) kaHist->Fill(p,dedx,1);
- if (w[2]>w[3] && w[2]>w[4] && w[2]>w[0]) piHist->Fill(p,dedx,1);
- if (w[0]>w[3] && w[0]>w[2] && w[0]>w[4]) elHist->Fill(p,dedx,1);
- }
-
- }
- delete event;
+ AliESDtrack *track=event->GetTrack(ntrk);
+ UInt_t status=track->GetStatus();
+
+ //select only tracks with the "combined PID"
+ if ((status&AliESDtrack::kESDpid)==0) continue;
+
+ Double_t w[10]; track->GetESDpid(w);
+ //count only "Kaon-like" tracks
+ if (w[3]>w[4] && w[3]>w[2] && w[3]>w[1] && w[3]>w[0]) nk++;
+ }
+ cout<<"Number of \"Kaon-like\" tracks : "<<nk<<endl;
+
+ //****** The loop over V0s
+ while (nv0--) {
+ AliESDv0 *v0=event->GetV0(nv0);
+ v0->ChangeMassHypothesis(310); // K0s
+ Double_t mass=v0->GetEffMass();
+ hm->Fill(mass);
+
+ Int_t pidx=v0->GetPindex(); // now let's get an access
+ AliESDtrack *track=event->GetTrack(pidx); // to the positive daughter
+ Double_t p=track->GetP();
+ hp->Fill(p);
+ }
+
+ //****** The loop over cascades
+ while (ncas--) {
+ AliESDcascade *cas=event->GetCascade(ncas);
+ Double_t q; //"quality" of the associated Lambda
+ cas->ChangeMassHypothesis(q,3312); // Xi-
+ // Here you do something with your Xis
+ // ...
+ // You can get the access to the daughters
+ }
+
}
+
+ delete event;
+ ef->Close();
+
timer.Stop(); timer.Print();
TCanvas *c1=new TCanvas("c1","",0,0,600,1200);
c1->Divide(1,2);
c1->cd(1);
- tpcHist->Draw();
- c1->cd(2);
- prHist->Draw();
- kaHist->Draw("same");
- piHist->Draw("same");
- elHist->Draw("same");
+ hm->Fit("gaus","","",V0mass-V0width,V0mass+V0width);
- ef->Close();
+ c1->cd(2);
+ hp->Fit("expo","","",0.3,2);
- return rc;
+ return 0;
}