]>
Commit | Line | Data |
---|---|---|
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 | ||
29 | extern AliRun *gAlice; | |
30 | extern TBenchmark *gBenchmark; | |
31 | extern TROOT *gROOT; | |
32 | ||
33 | static Int_t allpisel=0; | |
34 | static Int_t allkasel=0; | |
35 | static Int_t allprsel=0; | |
36 | static Int_t allnsel=0; | |
37 | ||
38 | ||
39 | Double_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 | ||
50 | Int_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 | } |