Updated version (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
c630aafd 8#if !defined( __CINT__) || defined(__MAKECINT__)
ae982df3 9 #include <Riostream.h>
10 #include "TKey.h"
11 #include "TFile.h"
c630aafd 12 #include "TH1F.h"
ae982df3 13 #include "TH2F.h"
14 #include "TCanvas.h"
15 #include "TStopwatch.h"
c630aafd 16 #include "TParticle.h"
17
18 #include "AliRun.h"
19 #include "AliRunLoader.h"
20 #include "AliLoader.h"
ae982df3 21
22 #include "AliESD.h"
23#endif
24
c630aafd 25extern AliRun *gAlice;
26
ae982df3 27Int_t AliESDanalysis(Int_t nev=1) {
8c6a71ab 28 TH2F *tpcHist=
c630aafd 29 new TH2F("tpcHist","TPC dE/dX vs momentum",100,0.,3.,100,0.,500.);
30 tpcHist->SetXTitle("p (GeV/c)"); tpcHist->SetYTitle("dE/dx (Arb. Units)");
8c6a71ab 31 tpcHist->SetMarkerStyle(8);
32 tpcHist->SetMarkerSize(0.3);
33
c630aafd 34 TH1F *piG=new TH1F("piG","",20,0.,3.);
35 TH1F *piR=new TH1F("piR","",20,0.,3.);
36 TH1F *piF=new TH1F("piF","",20,0.,3.);
37 TH1F *piGood=new TH1F("piGood","Combined PID for pions",20,0.,3.);
38 piGood->SetLineColor(4); piGood->SetXTitle("p (GeV/c)");
39 TH1F *piFake=new TH1F("piFake","",20,0.,3.); piFake->SetLineColor(2);
40
41 TH1F *kaG=new TH1F("kaG","",20,0.,3.);
42 TH1F *kaR=new TH1F("kaR","",20,0.,3.);
43 TH1F *kaF=new TH1F("kaF","",20,0.,3.);
44 TH1F *kaGood=new TH1F("kaGood","Combined PID for kaons",20,0.,3.);
45 kaGood->SetLineColor(4); kaGood->SetXTitle("p (GeV/c)");
46 TH1F *kaFake=new TH1F("kaFake","",20,0.,3.); kaFake->SetLineColor(2);
47
48 TH1F *prG=new TH1F("prG","",20,0.,3.);
49 TH1F *prR=new TH1F("prR","",20,0.,3.);
50 TH1F *prF=new TH1F("prF","",20,0.,3.);
51 TH1F *prGood=new TH1F("prGood","Combined PID for protons",20,0.,3.);
52 prGood->SetLineColor(4); prGood->SetXTitle("p (GeV/c)");
53 TH1F *prFake=new TH1F("prFake","",20,0.,3.); prFake->SetLineColor(2);
54
55 if (gAlice) {
56 delete gAlice->GetRunLoader();
57 delete gAlice;
58 gAlice=0;
59 }
60 AliRunLoader *rl = AliRunLoader::Open("galice.root");
61 if (rl == 0x0) {
62 cerr<<"Can not open session"<<endl;
63 return 1;
64 }
65 if (rl->LoadgAlice()) {
66 cerr<<"LoadgAlice returned error"<<endl;
67 delete rl;
68 return 1;
69 }
70 if (rl->LoadHeader()) {
71 cerr<<"LoadHeader returned error"<<endl;
72 delete rl;
73 return 1;
74 }
75 rl->LoadKinematics();
76
ae982df3 77 TFile *ef=TFile::Open("AliESDs.root");
78 if (!ef->IsOpen()) {cerr<<"Can't AliESDs.root !\n"; return 1;}
79
80 TStopwatch timer;
81 Int_t rc=0,n=0;
82 TKey *key=0;
83 TIter next(ef->GetListOfKeys());
84
8c6a71ab 85 //****** Tentative particle type "concentrations"
c630aafd 86 Double_t c[5]={0.01, 0.01, 0.85, 0.10, 0.05};
8c6a71ab 87
ae982df3 88 //******* The loop over events
89 while ((key=(TKey*)next())!=0) {
c630aafd 90 rl->GetEvent(n);
91
ae982df3 92 cerr<<"Processing event number : "<<n++<<endl;
c630aafd 93
ae982df3 94 AliESD *event=(AliESD*)key->ReadObj();
95
96 Int_t ntrk=event->GetNumberOfTracks();
97 cerr<<"Number of ESD tracks : "<<ntrk<<endl;
98 //****** The loop over tracks
c630aafd 99
100 Int_t pisel=0,kasel=0,prsel=0,nsel=0;
101
ae982df3 102 while (ntrk--) {
103 AliESDtrack *t=event->GetTrack(ntrk);
8c6a71ab 104
ae982df3 105 Double_t p=t->GetP();
8c6a71ab 106
ae982df3 107 if (t->GetStatus()&AliESDtrack::kTPCin) {
108 Double_t dedx=t->GetTPCsignal();
109 tpcHist->Fill(p,dedx,1);
110 }
8c6a71ab 111
c630aafd 112 UInt_t status=AliESDtrack::kESDpid;
113 status|=AliESDtrack::kTPCpid;
114 status|=AliESDtrack::kTOFpid;
115
116 if ((t->GetStatus()&status) == status) {
117 nsel++;
118
119 Int_t lab=TMath::Abs(t->GetLabel());
120 TParticle *part=gAlice->Particle(lab);
121 Int_t code=part->GetPdgCode();
122
123 Double_t r[10]; t->GetESDpid(r);
124
125 Double_t rcc=0.;
8c6a71ab 126 Int_t i;
c630aafd 127 for (i=0; i<AliESDtrack::kSPECIES; i++) rcc+=(c[i]*r[i]);
128 if (rcc==0.) continue;
8c6a71ab 129
130 //Here we apply Bayes' formula
131 Double_t w[10];
c630aafd 132 for (i=0; i<AliESDtrack::kSPECIES; i++) w[i]=c[i]*r[i]/rcc;
133
134 if (w[4]>w[3] && w[4]>w[2] && w[4]>w[1] && w[4]>w[0]) {//proton
135 prsel++;
136 prG->Fill(p);
137 if (TMath::Abs(code)==2212) prR->Fill(p);
138 else prF->Fill(p);
139 }
140
141 if (w[3]>w[4] && w[3]>w[2] && w[3]>w[1] && w[3]>w[0]) {//kaon
142 kasel++;
143 kaG->Fill(p);
144 if (TMath::Abs(code)==321) kaR->Fill(p);
145 else kaF->Fill(p);
146 }
147
148 if (w[2]>w[3] && w[2]>w[4] && w[2]>w[1] && w[2]>w[0]) {//pion
149 pisel++;
150 piG->Fill(p);
151 if (TMath::Abs(code)==211) piR->Fill(p);
152 else piF->Fill(p);
153 }
8c6a71ab 154
ae982df3 155 }
8c6a71ab 156
ae982df3 157 }
158 delete event;
c630aafd 159 cerr<<"Number of selected ESD tracks : "<<nsel<<endl;
160 cerr<<"Number of selected pion ESD tracks : "<<pisel<<endl;
161 cerr<<"Number of selected kaon ESD tracks : "<<kasel<<endl;
162 cerr<<"Number of selected proton ESD tracks : "<<prsel<<endl;
ae982df3 163 }
c630aafd 164
ae982df3 165 timer.Stop(); timer.Print();
166
167 TCanvas *c1=new TCanvas("c1","",0,0,600,1200);
c630aafd 168 c1->Divide(1,4);
ae982df3 169
170 c1->cd(1);
171 tpcHist->Draw();
c630aafd 172
ae982df3 173 c1->cd(2);
c630aafd 174 piG->Sumw2(); piF->Sumw2(); piR->Sumw2();
175 piGood->Divide(piR,piG,1,1,"b");
176 piFake->Divide(piF,piG,1,1,"b");
177 piGood->Draw("hist");
178 piFake->Draw("same");
179
180 c1->cd(3);
181 kaG->Sumw2(); kaF->Sumw2(); kaR->Sumw2();
182 kaGood->Divide(kaR,kaG,1,1,"b");
183 kaFake->Divide(kaF,kaG,1,1,"b");
184 kaGood->Draw("hist");
185 kaFake->Draw("same");
186
187 c1->cd(4);
188 prG->Sumw2(); prF->Sumw2(); prR->Sumw2();
189 prGood->Divide(prR,prG,1,1,"b");
190 prFake->Divide(prF,prG,1,1,"b");
191 prGood->Draw("hist");
192 prFake->Draw("same");
ae982df3 193
194 ef->Close();
195
c630aafd 196 delete rl;
197
ae982df3 198 return rc;
199}