1 //*************************************************************************
2 // Example of using the PID classes
3 // which calculates the efficiency and contamination of the combined PID.
5 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
6 //*************************************************************************
8 #if !defined( __CINT__) || defined(__MAKECINT__)
11 #include <Riostream.h>
15 #include <TParticle.h>
17 #include <TBenchmark.h>
23 #include <AliRunLoader.h>
25 #include <AliTPCpidESD.h>
26 #include <AliESDEvent.h>
29 extern AliRun *gAlice;
30 extern TBenchmark *gBenchmark;
33 static Int_t allpisel=0;
34 static Int_t allkasel=0;
35 static Int_t allprsel=0;
36 static Int_t allnsel=0;
39 Double_t tpcBethe(Double_t *x, Double_t *par) {
41 // This is the TPC Bethe-Bloch function normalised to 1 at the minimum
43 Double_t p=x[0]; //momentum
44 Double_t m=par[0]; //mass
46 return AliTPCpidESD::Bethe(bg);
50 Int_t AliPIDComparison(const Char_t *dir=".") {
51 gBenchmark->Start("AliPIDComparison");
53 ::Info("AliPIDComparison.C","Doing comparison...");
57 TH2F *tpcHist=(TH2F*)gROOT->FindObject("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);
64 const Char_t *hname[]={
65 "piG","piR","piF","piGood","piFake",
66 "kaG","kaR","kaF","kaGood","kaFake",
67 "prG","prR","prF","prGood","prFake"
69 Int_t nh=sizeof(hname)/sizeof(const Char_t *);
70 TH1F **hprt=new TH1F*[nh];
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();}
80 piGood->SetTitle("Combined PID for pions");
81 piGood->SetLineColor(4); piGood->SetXTitle("p (GeV/c)");
83 piFake->SetLineColor(2);
88 kaGood->SetTitle("Combined PID for kaons");
89 kaGood->SetLineColor(4); kaGood->SetXTitle("p (GeV/c)");
91 kaFake->SetLineColor(2);
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);
106 delete AliRunLoader::Instance();
110 sprintf(fname,"%s/galice.root",dir);
111 AliRunLoader *rl = AliRunLoader::Open(fname);
113 ::Error("AliPIDComparison.C","Can not open session !");
116 if (rl->LoadgAlice()) {
117 ::Error("AliPIDComparison.C","LoadgAlice returned error !");
121 if (rl->LoadHeader()) {
122 ::Error("AliPIDComparison.C","LoadHeader returned error !");
126 rl->LoadKinematics();
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 !");
135 AliESDEvent* event = new AliESDEvent();
136 TTree* tree = (TTree*) ef->Get("esdTree");
138 ::Error("AliPIDComparison.C", "no ESD tree found");
142 event->ReadFromTree(tree);
145 //****** Tentative particle type "concentrations"
146 Double_t c[5]={0.01, 0.01, 0.85, 0.10, 0.05};
147 AliPID::SetPriors(c);
150 //******* The loop over events
152 while (tree->GetEvent(e)) {
153 cout<<endl<<endl<<"********* Processing event number: "<<e<<"*******\n";
159 Int_t ntrk=event->GetNumberOfTracks();
160 cerr<<"Number of ESD tracks : "<<ntrk<<endl;
162 Int_t pisel=0,kasel=0,prsel=0,nsel=0;
164 AliStack *stack = rl->Stack();
167 AliESDtrack *t=event->GetTrack(ntrk);
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;
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;
181 Double_t p=t->GetInnerParam()->GetP();
182 Double_t dedx=t->GetTPCsignal()/47.;
183 tpcHist->Fill(p,dedx,1);
185 Int_t lab=TMath::Abs(t->GetLabel());
186 TParticle *part=stack->Particle(lab);
187 Int_t code=part->GetPdgCode();
189 Double_t r[10]; t->GetESDpid(r);
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);
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
204 if (TMath::Abs(code)!=2212) prF->Fill(p);
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
211 if (TMath::Abs(code)!=321) kaF->Fill(p);
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
218 if (TMath::Abs(code)!=211) piF->Fill(p);
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
225 if (TMath::Abs(code)!=11) piF->Fill(p);
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;
235 allnsel+=nsel; allpisel+=pisel; allkasel+=kasel; allprsel+=prsel;
237 } // ***** End of the loop over events
243 TCanvas *c1=(TCanvas*)gROOT->FindObject("c1");
245 c1=new TCanvas("c1","",0,0,600,1000);
252 TF1 *fpi=(TF1*)gROOT->FindObject("piBethe");
254 fpi=new TF1("piBethe",tpcBethe,pi,pa,1);
255 fpi->SetLineColor(2);
256 fpi->SetParameter(0,AliPID::ParticleMass(AliPID::kPion));
260 TF1 *fka=(TF1*)gROOT->FindObject("kaBethe");
262 fka=new TF1("kaBethe",tpcBethe,pi,pa,1);
263 fka->SetLineColor(3);
264 fka->SetParameter(0,AliPID::ParticleMass(AliPID::kKaon));
268 TF1 *fpr=(TF1*)gROOT->FindObject("prBethe");
270 fpr=new TF1("prBethe",tpcBethe,pi,pa,1);
271 fpr->SetLineColor(4);
272 fpr->SetParameter(0,AliPID::ParticleMass(AliPID::kProton));
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");
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");
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");
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;
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;
328 TFile fc("AliPIDComparison.root","RECREATE");
332 gBenchmark->Stop("AliPIDComparison");
333 gBenchmark->Show("AliPIDComparison");