First version of combined PID (Yu.Belikov)
[u/mrichter/AliRoot.git] / STEER / AliESDanalysis.C
CommitLineData
ae982df3 1//********************************************************************
8c6a71ab 2// Example (very naive for the moment) of the data analysis
ae982df3 3// using the ESD classes
8c6a71ab 4// It demonstrates the idea of the "combined PID".
ae982df3 5// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
6//********************************************************************
7
8#ifndef __CINT__
9 #include <Riostream.h>
10 #include "TKey.h"
11 #include "TFile.h"
12 #include "TH2F.h"
13 #include "TCanvas.h"
14 #include "TStopwatch.h"
15
16 #include "AliESD.h"
17#endif
18
19Int_t AliESDanalysis(Int_t nev=1) {
8c6a71ab 20 TH2F *tpcHist=
21 new TH2F("tpcHist","TPC dE/dX vs momentum",100,0.,4.,100,0.,500.);
22 tpcHist->SetXTitle("p (GeV/c)"); tpcHist->SetYTitle("dE/dx (Arb. Units)");
23 tpcHist->SetMarkerStyle(8);
24 tpcHist->SetMarkerSize(0.3);
25
26 TH2F *elHist=new TH2F("elHist","dE/dX vs momentum",100,0.,4.,50,0.,500.);
27 elHist->SetMarkerStyle(8);
28 elHist->SetMarkerSize(0.3);
29
30 TH2F *piHist=new TH2F("piHist","dE/dX vs momentum",100,0.,4.,50,0.,500.);
31 piHist->SetMarkerColor(2);
32 piHist->SetMarkerStyle(8);
33 piHist->SetMarkerSize(0.3);
34
35 TH2F *kaHist=new TH2F("kaHist","dE/dX vs momentum",100,0.,4.,100,0.,500.);
36 kaHist->SetMarkerColor(4);
37 kaHist->SetMarkerStyle(8);
38 kaHist->SetMarkerSize(0.3);
39
40 TH2F *prHist=new
41 TH2F("prHist","Classification into e, pi, K and p",100,0.,4.,100,0.,500.);
42 prHist->SetXTitle("p (GeV/c)"); prHist->SetYTitle("dE/dx (Arb. Units)");
43 prHist->SetMarkerColor(6);
44 prHist->SetMarkerStyle(8);
45 prHist->SetMarkerSize(0.3);
ae982df3 46
47 TFile *ef=TFile::Open("AliESDs.root");
48 if (!ef->IsOpen()) {cerr<<"Can't AliESDs.root !\n"; return 1;}
49
50 TStopwatch timer;
51 Int_t rc=0,n=0;
52 TKey *key=0;
53 TIter next(ef->GetListOfKeys());
54
8c6a71ab 55 //****** Tentative particle type "concentrations"
56 Double_t c[5]={0.05, 0., 0.85, 0.10, 0.05};
57
ae982df3 58 //******* The loop over events
59 while ((key=(TKey*)next())!=0) {
60 cerr<<"Processing event number : "<<n++<<endl;
61 AliESD *event=(AliESD*)key->ReadObj();
62
63 Int_t ntrk=event->GetNumberOfTracks();
64 cerr<<"Number of ESD tracks : "<<ntrk<<endl;
65 //****** The loop over tracks
66 while (ntrk--) {
67 AliESDtrack *t=event->GetTrack(ntrk);
8c6a71ab 68
ae982df3 69 Double_t p=t->GetP();
8c6a71ab 70
ae982df3 71 if (t->GetStatus()&AliESDtrack::kTPCin) {
72 Double_t dedx=t->GetTPCsignal();
73 tpcHist->Fill(p,dedx,1);
74 }
8c6a71ab 75
76 if (t->GetStatus()&AliESDtrack::kESDpid) {
77 Double_t dedx=t->GetTPCsignal();
78 Double_t r[10]; t->GetESDpid(r);
79 Double_t rc=0.;
80 Int_t i;
81 for (i=0; i<AliESDtrack::kSPECIES; i++) rc+=(c[i]*r[i]);
82 if (rc==0.) continue;
83
84 //Here we apply Bayes' formula
85 Double_t w[10];
86 for (i=0; i<AliESDtrack::kSPECIES; i++) w[i]=c[i]*r[i]/rc;
87
88 if (w[4]>w[3] && w[4]>w[2] && w[4]>w[0]) prHist->Fill(p,dedx,1);
89 if (w[3]>w[4] && w[3]>w[2] && w[3]>w[0]) kaHist->Fill(p,dedx,1);
90 if (w[2]>w[3] && w[2]>w[4] && w[2]>w[0]) piHist->Fill(p,dedx,1);
91 if (w[0]>w[3] && w[0]>w[2] && w[0]>w[4]) elHist->Fill(p,dedx,1);
ae982df3 92 }
8c6a71ab 93
ae982df3 94 }
95 delete event;
96 }
97 timer.Stop(); timer.Print();
98
99 TCanvas *c1=new TCanvas("c1","",0,0,600,1200);
100 c1->Divide(1,2);
101
102 c1->cd(1);
103 tpcHist->Draw();
104 c1->cd(2);
8c6a71ab 105 prHist->Draw();
106 kaHist->Draw("same");
107 piHist->Draw("same");
108 elHist->Draw("same");
ae982df3 109
110 ef->Close();
111
112 return rc;
113}