Bug corrected.
[u/mrichter/AliRoot.git] / STEER / AliPIDComparison.C
CommitLineData
8781a4f1 1//*************************************************************************
2// Example of using the PID classes
3// which calculates the efficiency and contamination of the combined PID.
4//
5// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
6//*************************************************************************
7
8#if !defined( __CINT__) || defined(__MAKECINT__)
9 #include <TMath.h>
10 #include <TError.h>
11 #include <Riostream.h>
12 #include <TF1.h>
13 #include <TH1F.h>
14 #include <TH2F.h>
15 #include <TParticle.h>
16 #include <TCanvas.h>
17 #include <TBenchmark.h>
18 #include <TFile.h>
19 #include <TTree.h>
20 #include <TROOT.h>
21
22 #include <AliStack.h>
23 #include <AliRunLoader.h>
24 #include <AliRun.h>
25 #include <AliTPCpidESD.h>
26 #include <AliESDEvent.h>
27#endif
28
29extern AliRun *gAlice;
30extern TBenchmark *gBenchmark;
31extern TROOT *gROOT;
32
33static Int_t allpisel=0;
34static Int_t allkasel=0;
35static Int_t allprsel=0;
36static Int_t allnsel=0;
37
38
39Double_t tpcBethe(Double_t *x, Double_t *par) {
40 //
41 // This is the TPC Bethe-Bloch function normalised to 1 at the minimum
42 //
43 Double_t p=x[0]; //momentum
44 Double_t m=par[0]; //mass
45 Double_t bg=p/m;
46 return AliTPCpidESD::Bethe(bg);
47}
48
49
50Int_t AliPIDComparison(const Char_t *dir=".") {
51 gBenchmark->Start("AliPIDComparison");
52
53 ::Info("AliPIDComparison.C","Doing comparison...");
54
55 Double_t pi=0.2,pa=3;
56
57 TH2F *tpcHist=(TH2F*)gROOT->FindObject("tpcHist");
58 if (!tpcHist)
59 tpcHist=new TH2F("tpcHist","TPC dE/dX vs momentum",100,pi,pa,100,0.,6.);
60 tpcHist->SetXTitle("p (GeV/c)"); tpcHist->SetYTitle("dE/dx (Arb. Units)");
61 tpcHist->SetMarkerStyle(8);
62 tpcHist->SetMarkerSize(0.3);
63
64 const Char_t *hname[]={
65 "piG","piR","piF","piGood","piFake",
66 "kaG","kaR","kaF","kaGood","kaFake",
67 "prG","prR","prF","prGood","prFake"
68 };
69 Int_t nh=sizeof(hname)/sizeof(const Char_t *);
70 TH1F **hprt=new TH1F*[nh];
71
72 for (Int_t i=0; i<nh; i++) {
73 hprt[i]=(TH1F*)gROOT->FindObject(hname[i]);
74 if (hprt[i]==0) {hprt[i]=new TH1F(hname[i],"",20,pi,pa);hprt[i]->Sumw2();}
75 }
76 TH1F *piG=hprt[0];
77 TH1F *piR=hprt[1];
78 TH1F *piF=hprt[2];
79 TH1F *piGood=hprt[3];
80 piGood->SetTitle("Combined PID for pions");
81 piGood->SetLineColor(4); piGood->SetXTitle("p (GeV/c)");
82 TH1F *piFake=hprt[4];
83 piFake->SetLineColor(2);
84 TH1F *kaG=hprt[5];
85 TH1F *kaR=hprt[6];
86 TH1F *kaF=hprt[7];
87 TH1F *kaGood=hprt[8];
88 kaGood->SetTitle("Combined PID for kaons");
89 kaGood->SetLineColor(4); kaGood->SetXTitle("p (GeV/c)");
90 TH1F *kaFake=hprt[9];
91 kaFake->SetLineColor(2);
92 TH1F *prG=hprt[10];
93 TH1F *prR=hprt[11];
94 TH1F *prF=hprt[12];
95 TH1F *prGood=hprt[13];
96 prGood->SetTitle("Combined PID for protons");
97 prGood->SetLineColor(4); prGood->SetXTitle("p (GeV/c)");
98 TH1F *prFake=hprt[14];
99 prFake->SetLineColor(2);
100
101 delete[] hprt;
102
103 Char_t fname[100];
104
105 if (gAlice) {
33c3c91a 106 delete AliRunLoader::Instance();
8781a4f1 107 delete gAlice;
108 gAlice=0;
109 }
110 sprintf(fname,"%s/galice.root",dir);
111 AliRunLoader *rl = AliRunLoader::Open(fname);
112 if (rl == 0x0) {
113 ::Error("AliPIDComparison.C","Can not open session !");
114 return 1;
115 }
116 if (rl->LoadgAlice()) {
117 ::Error("AliPIDComparison.C","LoadgAlice returned error !");
118 delete rl;
119 return 1;
120 }
121 if (rl->LoadHeader()) {
122 ::Error("AliPIDComparison.C","LoadHeader returned error !");
123 delete rl;
124 return 1;
125 }
126 rl->LoadKinematics();
127
128 sprintf(fname,"%s/AliESDs.root",dir);
129 TFile *ef=TFile::Open(fname);
130 if (!ef || !ef->IsOpen()) {
131 ::Error("AliPIDComparison.C","Can't AliESDs.root !");
132 delete rl;
133 return 1;
134 }
135 AliESDEvent* event = new AliESDEvent();
136 TTree* tree = (TTree*) ef->Get("esdTree");
137 if (!tree) {
138 ::Error("AliPIDComparison.C", "no ESD tree found");
139 delete rl;
140 return 1;
141 }
142 event->ReadFromTree(tree);
143
144
145 //****** Tentative particle type "concentrations"
146 Double_t c[5]={0.01, 0.01, 0.85, 0.10, 0.05};
147 AliPID::SetPriors(c);
148
149
150 //******* The loop over events
151 Int_t e=0;
152 while (tree->GetEvent(e)) {
153 cout<<endl<<endl<<"********* Processing event number: "<<e<<"*******\n";
154
155 rl->GetEvent(e);
156
157 e++;
158
159 Int_t ntrk=event->GetNumberOfTracks();
160 cerr<<"Number of ESD tracks : "<<ntrk<<endl;
161
162 Int_t pisel=0,kasel=0,prsel=0,nsel=0;
163
164 AliStack *stack = rl->Stack();
165
166 while (ntrk--) {
167 AliESDtrack *t=event->GetTrack(ntrk);
168
169 //*** Some track quality cuts ****
170 if (t->GetITSclusters(0) < 6 ) continue;
171 if (t->GetTPCclusters(0) < 60) continue;
172 //if (t->GetTRDclusters(0) < 60) continue;
173
174 if (!t->IsOn(AliESDtrack::kITSpid)) continue;
175 if (!t->IsOn(AliESDtrack::kTPCpid)) continue;
176 //if (!t->IsOn(AliESDtrack::kTRDpid)) continue;
177 if (!t->IsOn(AliESDtrack::kTOFpid)) continue;
178 {
179 nsel++;
180
181 Double_t p=t->GetInnerParam()->GetP();
182 Double_t dedx=t->GetTPCsignal()/47.;
183 tpcHist->Fill(p,dedx,1);
184
185 Int_t lab=TMath::Abs(t->GetLabel());
186 TParticle *part=stack->Particle(lab);
187 Int_t code=part->GetPdgCode();
188
189 Double_t r[10]; t->GetESDpid(r);
190
191 AliPID pid(r);
192
193 Double_t w[10];
194 w[0]=pid.GetProbability(AliPID::kElectron);
195 w[1]=pid.GetProbability(AliPID::kMuon);
196 w[2]=pid.GetProbability(AliPID::kPion);
197 w[3]=pid.GetProbability(AliPID::kKaon);
198 w[4]=pid.GetProbability(AliPID::kProton);
199
200 if (TMath::Abs(code)==2212) prR->Fill(p);
201 if (w[4]>w[3] && w[4]>w[2] && w[4]>w[1] && w[4]>w[0]) {//proton
202 prsel++;
203 prG->Fill(p);
204 if (TMath::Abs(code)!=2212) prF->Fill(p);
205 }
206
207 if (TMath::Abs(code)==321) kaR->Fill(p);
208 if (w[3]>w[4] && w[3]>w[2] && w[3]>w[1] && w[3]>w[0]) {//kaon
209 kasel++;
210 kaG->Fill(p);
211 if (TMath::Abs(code)!=321) kaF->Fill(p);
212 }
213
214 if (TMath::Abs(code)==211) piR->Fill(p);
215 if (w[2]>w[4] && w[2]>w[3] && w[2]>w[0] && w[2]>w[1]) {//pion
216 pisel++;
217 piG->Fill(p);
218 if (TMath::Abs(code)!=211) piF->Fill(p);
219 }
220 /*
221 if (TMath::Abs(code)==11) piR->Fill(p);
222 if (w[0]>w[4] && w[0]>w[3] && w[0]>w[3] && w[0]>w[1]) {//electron
223 pisel++;
224 piG->Fill(p);
225 if (TMath::Abs(code)!=11) piF->Fill(p);
226 }
227 */
228 }
229 }
230 cout<<"Number of selected ESD tracks : "<<nsel<<endl;
231 cout<<"Number of selected pion ESD tracks : "<<pisel<<endl;
232 cout<<"Number of selected kaon ESD tracks : "<<kasel<<endl;
233 cout<<"Number of selected proton ESD tracks : "<<prsel<<endl;
234
235 allnsel+=nsel; allpisel+=pisel; allkasel+=kasel; allprsel+=prsel;
236
237 } // ***** End of the loop over events
238
239 delete event;
240 delete tree;
241 ef->Close();
242
243 TCanvas *c1=(TCanvas*)gROOT->FindObject("c1");
244 if (!c1) {
245 c1=new TCanvas("c1","",0,0,600,1000);
246 c1->Divide(1,4);
247 }
248
249 c1->cd(1);
250 tpcHist->Draw();
251
252 TF1 *fpi=(TF1*)gROOT->FindObject("piBethe");
253 if (!fpi) {
254 fpi=new TF1("piBethe",tpcBethe,pi,pa,1);
255 fpi->SetLineColor(2);
256 fpi->SetParameter(0,AliPID::ParticleMass(AliPID::kPion));
257 }
258 fpi->Draw("same");
259
260 TF1 *fka=(TF1*)gROOT->FindObject("kaBethe");
261 if (!fka) {
262 fka=new TF1("kaBethe",tpcBethe,pi,pa,1);
263 fka->SetLineColor(3);
264 fka->SetParameter(0,AliPID::ParticleMass(AliPID::kKaon));
265 }
266 fka->Draw("same");
267
268 TF1 *fpr=(TF1*)gROOT->FindObject("prBethe");
269 if (!fpr) {
270 fpr=new TF1("prBethe",tpcBethe,pi,pa,1);
271 fpr->SetLineColor(4);
272 fpr->SetParameter(0,AliPID::ParticleMass(AliPID::kProton));
273 }
274 fpr->Draw("same");
275
276
277 c1->cd(2);
278 //piG->Sumw2(); piF->Sumw2(); piR->Sumw2();
279 piGood->Add(piG,piF,1,-1);
280 piGood->Divide(piGood,piR,1,1,"b");
281 piFake->Divide(piF,piG,1,1,"b");
282 piGood->SetMinimum(0);
283 piGood->Draw("hist");
284 piFake->Draw("same");
285
286 c1->cd(3);
287 //kaG->Sumw2(); kaF->Sumw2(); kaR->Sumw2();
288 kaGood->Add(kaG,kaF,1,-1);
289 kaGood->Divide(kaGood,kaR,1,1,"b");
290 kaFake->Divide(kaF,kaG,1,1,"b");
291 kaGood->SetMinimum(0);
292 kaGood->Draw("hist");
293 kaFake->Draw("same");
294
295 c1->cd(4);
296 //prG->Sumw2(); prF->Sumw2(); prR->Sumw2();
297 prGood->Add(prG,prF,1,-1);
298 prGood->Divide(prGood,prR,1,1,"b");
299 prFake->Divide(prF,prG,1,1,"b");
300 prGood->SetMinimum(0);
301 prGood->Draw("hist");
302 prFake->Draw("same");
303
304 c1->Update();
305
306 cout<<endl;
307 cout<<endl;
308 e=(Int_t)piR->GetEntries();
309 Int_t o=(Int_t)piG->GetEntries();
310 if (e*o) cout<<"Efficiency (contamination) for pions : "<<
311 piGood->GetEntries()/e<<'('<<piF->GetEntries()/o<<')'<<endl;
312 e=(Int_t)kaR->GetEntries();
313 o=(Int_t)kaG->GetEntries();
314 if (e*o) cout<<"Efficiency (contamination) for kaons : "<<
315 kaGood->GetEntries()/e<<'('<<kaF->GetEntries()/o<<')'<<endl;
316 e=(Int_t)prR->GetEntries();
317 o=(Int_t)prG->GetEntries();
318 if (e*o) cout<<"Efficiency (contamination) for protons : "<<
319 prGood->GetEntries()/e<<'('<<prF->GetEntries()/o<<')'<<endl;
320 cout<<endl;
321
322 cout<<"Total number of selected ESD tracks : "<<allnsel<<endl;
323 cout<<"Total number of selected pion ESD tracks : "<<allpisel<<endl;
324 cout<<"Total number of selected kaon ESD tracks : "<<allkasel<<endl;
325 cout<<"Total number of selected proton ESD tracks : "<<allprsel<<endl;
326
327 ef->Close();
328 TFile fc("AliPIDComparison.root","RECREATE");
329 c1->Write();
330 fc.Close();
331
332 gBenchmark->Stop("AliPIDComparison");
333 gBenchmark->Show("AliPIDComparison");
334
335 delete rl;
336 return 0;
337}