1 //********************************************************************
2 // Example (very naive for the moment) of using the ESD classes
3 // It demonstrates the idea of the "combined PID".
4 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
5 //********************************************************************
7 #if !defined( __CINT__) || defined(__MAKECINT__)
10 #include <Riostream.h>
13 #include <TParticle.h>
15 #include <TBenchmark.h>
21 #include <AliRunLoader.h>
23 #include <AliESDEvent.h>
26 extern AliRun *gAlice;
27 extern TBenchmark *gBenchmark;
30 static Int_t allpisel=0;
31 static Int_t allkasel=0;
32 static Int_t allprsel=0;
33 static Int_t allnsel=0;
35 Int_t AliESDComparison(const Char_t *dir=".") {
36 gBenchmark->Start("AliESDComparison");
38 ::Info("AliESDComparison.C","Doing comparison...");
42 TH2F *tpcHist=(TH2F*)gROOT->FindObject("tpcHist");
44 tpcHist=new TH2F("tpcHist","TPC dE/dX vs momentum",100,pi,pa,100,0.,300.);
45 tpcHist->SetXTitle("p (GeV/c)"); tpcHist->SetYTitle("dE/dx (Arb. Units)");
46 tpcHist->SetMarkerStyle(8);
47 tpcHist->SetMarkerSize(0.3);
49 const Char_t *hname[]={
50 "piG","piR","piF","piGood","piFake",
51 "kaG","kaR","kaF","kaGood","kaFake",
52 "prG","prR","prF","prGood","prFake"
54 Int_t nh=sizeof(hname)/sizeof(const Char_t *);
55 TH1F **hprt=new TH1F*[nh];
57 for (Int_t i=0; i<nh; i++) {
58 hprt[i]=(TH1F*)gROOT->FindObject(hname[i]);
59 if (hprt[i]==0) {hprt[i]=new TH1F(hname[i],"",20,pi,pa);hprt[i]->Sumw2();}
65 piGood->SetTitle("Combined PID for pions");
66 piGood->SetLineColor(4); piGood->SetXTitle("p (GeV/c)");
68 piFake->SetLineColor(2);
73 kaGood->SetTitle("Combined PID for kaons");
74 kaGood->SetLineColor(4); kaGood->SetXTitle("p (GeV/c)");
76 kaFake->SetLineColor(2);
80 TH1F *prGood=hprt[13];
81 prGood->SetTitle("Combined PID for protons");
82 prGood->SetLineColor(4); prGood->SetXTitle("p (GeV/c)");
83 TH1F *prFake=hprt[14];
84 prFake->SetLineColor(2);
91 delete gAlice->GetRunLoader();
95 sprintf(fname,"%s/galice.root",dir);
96 AliRunLoader *rl = AliRunLoader::Open(fname);
98 ::Error("AliESDComparison.C","Can not open session !");
101 if (rl->LoadgAlice()) {
102 ::Error("AliESDComparison.C","LoadgAlice returned error !");
106 if (rl->LoadHeader()) {
107 ::Error("AliESDComparison.C","LoadHeader returned error !");
111 rl->LoadKinematics();
113 sprintf(fname,"%s/AliESDs.root",dir);
114 TFile *ef=TFile::Open(fname);
115 if (!ef || !ef->IsOpen()) {
116 ::Error("AliESDComparison.C","Can't AliESDs.root !");
120 AliESDEvent* event = new AliESDEvent();
121 TTree* tree = (TTree*) ef->Get("esdTree");
123 ::Error("AliESDComparison.C", "no ESD tree found");
127 event->ReadFromTree(tree);
129 //****** Tentative particle type "concentrations"
130 Double_t c[5]={0.01, 0.01, 0.85, 0.10, 0.05};
131 //Double_t c[5]={0.2, 0.2, 0.2, 0.2, 0.2};
132 AliPID::SetPriors(c);
135 //******* The loop over events
137 while (tree->GetEvent(e)) {
138 cout<<endl<<endl<<"********* Processing event number: "<<e<<"*******\n";
144 Int_t ntrk=event->GetNumberOfTracks();
145 cerr<<"Number of ESD tracks : "<<ntrk<<endl;
147 Int_t pisel=0,kasel=0,prsel=0,nsel=0;
149 AliStack *stack = rl->Stack();
152 AliESDtrack *t=event->GetTrack(ntrk);
154 //*** Some track quality cuts ****
155 if (t->GetITSclusters(0) < 6 ) continue;
156 if (t->GetTPCclusters(0) < 60) continue;
157 //if (t->GetTRDclusters(0) < 60) continue;
159 if (!t->IsOn(AliESDtrack::kESDpid)) continue;
160 if (!t->IsOn(AliESDtrack::kITSpid)) continue;
161 if (!t->IsOn(AliESDtrack::kTPCpid)) continue;
162 //if (!t->IsOn(AliESDtrack::kTRDpid)) continue;
163 if (!t->IsOn(AliESDtrack::kTOFpid)) continue;
167 Double_t p=t->GetP();
168 Double_t dedx=t->GetTPCsignal();
169 tpcHist->Fill(p,dedx,1);
171 Int_t lab=TMath::Abs(t->GetLabel());
172 TParticle *part=stack->Particle(lab);
173 Int_t code=part->GetPdgCode();
175 Double_t r[10]; t->GetESDpid(r);
184 w[0]=pid.GetProbability(AliPID::kElectron);
185 w[1]=pid.GetProbability(AliPID::kMuon);
186 w[2]=pid.GetProbability(AliPID::kPion);
187 w[3]=pid.GetProbability(AliPID::kKaon);
188 w[4]=pid.GetProbability(AliPID::kProton);
191 if (TMath::Abs(code)==2212) prR->Fill(p);
192 if (w[4]>w[3] && w[4]>w[2] && w[4]>w[1] && w[4]>w[0]) {//proton
195 if (TMath::Abs(code)!=2212) prF->Fill(p);
198 if (TMath::Abs(code)==321) kaR->Fill(p);
199 if (w[3]>w[4] && w[3]>w[2] && w[3]>w[1] && w[3]>w[0]) {//kaon
202 if (TMath::Abs(code)!=321) kaF->Fill(p);
205 if (TMath::Abs(code)==211) piR->Fill(p);
206 if (w[2]>w[4] && w[2]>w[3] && w[2]>w[0] && w[2]>w[1]) {//pion
209 if (TMath::Abs(code)!=211) piF->Fill(p);
213 cout<<"Number of selected ESD tracks : "<<nsel<<endl;
214 cout<<"Number of selected pion ESD tracks : "<<pisel<<endl;
215 cout<<"Number of selected kaon ESD tracks : "<<kasel<<endl;
216 cout<<"Number of selected proton ESD tracks : "<<prsel<<endl;
218 allnsel+=nsel; allpisel+=pisel; allkasel+=kasel; allprsel+=prsel;
220 } // ***** End of the loop over events
226 TCanvas *c1=(TCanvas*)gROOT->FindObject("c1");
228 c1=new TCanvas("c1","",0,0,600,1200);
236 //piG->Sumw2(); piF->Sumw2(); piR->Sumw2();
237 piGood->Add(piG,piF,1,-1);
238 piGood->Divide(piGood,piR,1,1,"b");
239 piFake->Divide(piF,piG,1,1,"b");
240 piGood->Draw("hist");
241 piFake->Draw("same");
244 //kaG->Sumw2(); kaF->Sumw2(); kaR->Sumw2();
245 kaGood->Add(kaG,kaF,1,-1);
246 kaGood->Divide(kaGood,kaR,1,1,"b");
247 kaFake->Divide(kaF,kaG,1,1,"b");
248 kaGood->Draw("hist");
249 kaFake->Draw("same");
252 //prG->Sumw2(); prF->Sumw2(); prR->Sumw2();
253 prGood->Add(prG,prF,1,-1);
254 prGood->Divide(prGood,prR,1,1,"b");
255 prFake->Divide(prF,prG,1,1,"b");
256 prGood->Draw("hist");
257 prFake->Draw("same");
263 e=(Int_t)piR->GetEntries();
264 Int_t o=(Int_t)piG->GetEntries();
265 if (e*o) cout<<"Efficiency (contamination) for pions : "<<
266 piGood->GetEntries()/e<<'('<<piF->GetEntries()/o<<')'<<endl;
267 e=(Int_t)kaR->GetEntries();
268 o=(Int_t)kaG->GetEntries();
269 if (e*o) cout<<"Efficiency (contamination) for kaons : "<<
270 kaGood->GetEntries()/e<<'('<<kaF->GetEntries()/o<<')'<<endl;
271 e=(Int_t)prR->GetEntries();
272 o=(Int_t)prG->GetEntries();
273 if (e*o) cout<<"Efficiency (contamination) for protons : "<<
274 prGood->GetEntries()/e<<'('<<prF->GetEntries()/o<<')'<<endl;
277 cout<<"Total number of selected ESD tracks : "<<allnsel<<endl;
278 cout<<"Total number of selected pion ESD tracks : "<<allpisel<<endl;
279 cout<<"Total number of selected kaon ESD tracks : "<<allkasel<<endl;
280 cout<<"Total number of selected proton ESD tracks : "<<allprsel<<endl;
283 TFile fc("AliESDComparison.root","RECREATE");
287 gBenchmark->Stop("AliESDComparison");
288 gBenchmark->Show("AliESDComparison");