]>
Commit | Line | Data |
---|---|---|
0717295b | 1 | //******************************************************************** |
2 | // Example (very naive for the moment) of using the ESD classes | |
42345b19 | 3 | // It demonstrates the idea of the "combined PID". |
4 | // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch | |
0717295b | 5 | //******************************************************************** |
6 | ||
7 | #if !defined( __CINT__) || defined(__MAKECINT__) | |
4fa1e629 | 8 | #include <TMath.h> |
9 | #include <TError.h> | |
0717295b | 10 | #include <Riostream.h> |
4fa1e629 | 11 | #include <TH1F.h> |
12 | #include <TH2F.h> | |
13 | #include <TParticle.h> | |
14 | #include <TCanvas.h> | |
15 | #include <TBenchmark.h> | |
8b462fd8 | 16 | #include <TFile.h> |
17 | #include <TTree.h> | |
4fa1e629 | 18 | #include <TROOT.h> |
19 | ||
20 | #include <AliStack.h> | |
21 | #include <AliRunLoader.h> | |
22 | #include <AliRun.h> | |
86ad5fcb | 23 | #include <AliESDEvent.h> |
0717295b | 24 | #endif |
25 | ||
26 | extern AliRun *gAlice; | |
4fa1e629 | 27 | extern TBenchmark *gBenchmark; |
42345b19 | 28 | extern TROOT *gROOT; |
0717295b | 29 | |
4fa1e629 | 30 | static Int_t allpisel=0; |
31 | static Int_t allkasel=0; | |
32 | static Int_t allprsel=0; | |
33 | static Int_t allnsel=0; | |
34 | ||
3f7a17bb | 35 | Int_t AliESDComparison(const Char_t *dir=".") { |
4fa1e629 | 36 | gBenchmark->Start("AliESDComparison"); |
37 | ||
38 | ::Info("AliESDComparison.C","Doing comparison..."); | |
39 | ||
40 | Double_t pi=0.2,pa=3; | |
f2eae663 | 41 | |
42345b19 | 42 | TH2F *tpcHist=(TH2F*)gROOT->FindObject("tpcHist"); |
43 | if (!tpcHist) | |
f2eae663 | 44 | tpcHist=new TH2F("tpcHist","TPC dE/dX vs momentum",100,pi,pa,100,0.,300.); |
0717295b | 45 | tpcHist->SetXTitle("p (GeV/c)"); tpcHist->SetYTitle("dE/dx (Arb. Units)"); |
42345b19 | 46 | tpcHist->SetMarkerStyle(8); |
47 | tpcHist->SetMarkerSize(0.3); | |
0717295b | 48 | |
42345b19 | 49 | const Char_t *hname[]={ |
50 | "piG","piR","piF","piGood","piFake", | |
51 | "kaG","kaR","kaF","kaGood","kaFake", | |
52 | "prG","prR","prF","prGood","prFake" | |
53 | }; | |
54 | Int_t nh=sizeof(hname)/sizeof(const Char_t *); | |
55 | TH1F **hprt=new TH1F*[nh]; | |
56 | ||
57 | for (Int_t i=0; i<nh; i++) { | |
58 | hprt[i]=(TH1F*)gROOT->FindObject(hname[i]); | |
f2eae663 | 59 | if (hprt[i]==0) {hprt[i]=new TH1F(hname[i],"",20,pi,pa);hprt[i]->Sumw2();} |
42345b19 | 60 | } |
61 | TH1F *piG=hprt[0]; | |
62 | TH1F *piR=hprt[1]; | |
63 | TH1F *piF=hprt[2]; | |
64 | TH1F *piGood=hprt[3]; | |
65 | piGood->SetTitle("Combined PID for pions"); | |
66 | piGood->SetLineColor(4); piGood->SetXTitle("p (GeV/c)"); | |
67 | TH1F *piFake=hprt[4]; | |
68 | piFake->SetLineColor(2); | |
69 | TH1F *kaG=hprt[5]; | |
70 | TH1F *kaR=hprt[6]; | |
71 | TH1F *kaF=hprt[7]; | |
72 | TH1F *kaGood=hprt[8]; | |
73 | kaGood->SetTitle("Combined PID for kaons"); | |
74 | kaGood->SetLineColor(4); kaGood->SetXTitle("p (GeV/c)"); | |
75 | TH1F *kaFake=hprt[9]; | |
76 | kaFake->SetLineColor(2); | |
77 | TH1F *prG=hprt[10]; | |
78 | TH1F *prR=hprt[11]; | |
79 | TH1F *prF=hprt[12]; | |
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); | |
85 | ||
86 | delete[] hprt; | |
87 | ||
88 | Char_t fname[100]; | |
0717295b | 89 | |
90 | if (gAlice) { | |
91 | delete gAlice->GetRunLoader(); | |
92 | delete gAlice; | |
93 | gAlice=0; | |
94 | } | |
42345b19 | 95 | sprintf(fname,"%s/galice.root",dir); |
96 | AliRunLoader *rl = AliRunLoader::Open(fname); | |
0717295b | 97 | if (rl == 0x0) { |
4fa1e629 | 98 | ::Error("AliESDComparison.C","Can not open session !"); |
0717295b | 99 | return 1; |
100 | } | |
101 | if (rl->LoadgAlice()) { | |
4fa1e629 | 102 | ::Error("AliESDComparison.C","LoadgAlice returned error !"); |
0717295b | 103 | delete rl; |
104 | return 1; | |
105 | } | |
106 | if (rl->LoadHeader()) { | |
4fa1e629 | 107 | ::Error("AliESDComparison.C","LoadHeader returned error !"); |
0717295b | 108 | delete rl; |
109 | return 1; | |
110 | } | |
111 | rl->LoadKinematics(); | |
0717295b | 112 | |
42345b19 | 113 | sprintf(fname,"%s/AliESDs.root",dir); |
114 | TFile *ef=TFile::Open(fname); | |
4fa1e629 | 115 | if (!ef || !ef->IsOpen()) { |
116 | ::Error("AliESDComparison.C","Can't AliESDs.root !"); | |
117 | delete rl; | |
118 | return 1; | |
119 | } | |
86ad5fcb | 120 | AliESDEvent* event = new AliESDEvent(); |
8b462fd8 | 121 | TTree* tree = (TTree*) ef->Get("esdTree"); |
122 | if (!tree) { | |
123 | ::Error("AliESDComparison.C", "no ESD tree found"); | |
124 | delete rl; | |
125 | return 1; | |
126 | } | |
86ad5fcb | 127 | event->ReadFromTree(tree); |
0717295b | 128 | |
129 | //****** Tentative particle type "concentrations" | |
130 | Double_t c[5]={0.01, 0.01, 0.85, 0.10, 0.05}; | |
c2be14cf | 131 | //Double_t c[5]={0.2, 0.2, 0.2, 0.2, 0.2}; |
132 | AliPID::SetPriors(c); | |
133 | ||
0717295b | 134 | |
135 | //******* The loop over events | |
4fa1e629 | 136 | Int_t e=0; |
8b462fd8 | 137 | while (tree->GetEvent(e)) { |
4fa1e629 | 138 | cout<<endl<<endl<<"********* Processing event number: "<<e<<"*******\n"; |
0717295b | 139 | |
8b462fd8 | 140 | rl->GetEvent(e); |
4fa1e629 | 141 | |
142 | e++; | |
0717295b | 143 | |
4fa1e629 | 144 | Int_t ntrk=event->GetNumberOfTracks(); |
145 | cerr<<"Number of ESD tracks : "<<ntrk<<endl; | |
0717295b | 146 | |
4fa1e629 | 147 | Int_t pisel=0,kasel=0,prsel=0,nsel=0; |
0717295b | 148 | |
4fa1e629 | 149 | AliStack *stack = rl->Stack(); |
0717295b | 150 | |
4fa1e629 | 151 | while (ntrk--) { |
1efd642f | 152 | AliESDtrack *t=event->GetTrack(ntrk); |
153 | ||
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; | |
158 | ||
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; | |
164 | { | |
4fa1e629 | 165 | nsel++; |
f2eae663 | 166 | |
4fa1e629 | 167 | Double_t p=t->GetP(); |
168 | Double_t dedx=t->GetTPCsignal(); | |
169 | tpcHist->Fill(p,dedx,1); | |
0717295b | 170 | |
4fa1e629 | 171 | Int_t lab=TMath::Abs(t->GetLabel()); |
172 | TParticle *part=stack->Particle(lab); | |
173 | Int_t code=part->GetPdgCode(); | |
0717295b | 174 | |
4fa1e629 | 175 | Double_t r[10]; t->GetESDpid(r); |
1efd642f | 176 | //t->GetITSpid(r); |
c2be14cf | 177 | //t->GetTPCpid(r); |
1efd642f | 178 | //t->GetTRDpid(r); |
179 | //t->GetTOFpid(r); | |
0717295b | 180 | |
c2be14cf | 181 | AliPID pid(r); |
0717295b | 182 | |
4fa1e629 | 183 | Double_t w[10]; |
c2be14cf | 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); | |
189 | ||
0717295b | 190 | |
4fa1e629 | 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 | |
193 | prsel++; | |
194 | prG->Fill(p); | |
195 | if (TMath::Abs(code)!=2212) prF->Fill(p); | |
196 | } | |
0717295b | 197 | |
4fa1e629 | 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 | |
200 | kasel++; | |
201 | kaG->Fill(p); | |
202 | if (TMath::Abs(code)!=321) kaF->Fill(p); | |
203 | } | |
0717295b | 204 | |
4fa1e629 | 205 | if (TMath::Abs(code)==211) piR->Fill(p); |
c2be14cf | 206 | if (w[2]>w[4] && w[2]>w[3] && w[2]>w[0] && w[2]>w[1]) {//pion |
4fa1e629 | 207 | pisel++; |
208 | piG->Fill(p); | |
209 | if (TMath::Abs(code)!=211) piF->Fill(p); | |
210 | } | |
211 | } | |
212 | } | |
4fa1e629 | 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; | |
0717295b | 217 | |
4fa1e629 | 218 | allnsel+=nsel; allpisel+=pisel; allkasel+=kasel; allprsel+=prsel; |
f2eae663 | 219 | |
4fa1e629 | 220 | } // ***** End of the loop over events |
0717295b | 221 | |
8b462fd8 | 222 | delete event; |
4514157e | 223 | delete tree; |
4fa1e629 | 224 | ef->Close(); |
0717295b | 225 | |
42345b19 | 226 | TCanvas *c1=(TCanvas*)gROOT->FindObject("c1"); |
4fa1e629 | 227 | if (!c1) { |
228 | c1=new TCanvas("c1","",0,0,600,1200); | |
229 | c1->Divide(1,4); | |
230 | } | |
0717295b | 231 | |
232 | c1->cd(1); | |
233 | tpcHist->Draw(); | |
234 | ||
235 | c1->cd(2); | |
42345b19 | 236 | //piG->Sumw2(); piF->Sumw2(); piR->Sumw2(); |
f2eae663 | 237 | piGood->Add(piG,piF,1,-1); |
238 | piGood->Divide(piGood,piR,1,1,"b"); | |
0717295b | 239 | piFake->Divide(piF,piG,1,1,"b"); |
240 | piGood->Draw("hist"); | |
241 | piFake->Draw("same"); | |
242 | ||
243 | c1->cd(3); | |
42345b19 | 244 | //kaG->Sumw2(); kaF->Sumw2(); kaR->Sumw2(); |
f2eae663 | 245 | kaGood->Add(kaG,kaF,1,-1); |
246 | kaGood->Divide(kaGood,kaR,1,1,"b"); | |
0717295b | 247 | kaFake->Divide(kaF,kaG,1,1,"b"); |
248 | kaGood->Draw("hist"); | |
249 | kaFake->Draw("same"); | |
250 | ||
251 | c1->cd(4); | |
42345b19 | 252 | //prG->Sumw2(); prF->Sumw2(); prR->Sumw2(); |
f2eae663 | 253 | prGood->Add(prG,prF,1,-1); |
254 | prGood->Divide(prGood,prR,1,1,"b"); | |
0717295b | 255 | prFake->Divide(prF,prG,1,1,"b"); |
256 | prGood->Draw("hist"); | |
257 | prFake->Draw("same"); | |
258 | ||
4fa1e629 | 259 | c1->Update(); |
260 | ||
261 | cout<<endl; | |
262 | cout<<endl; | |
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; | |
275 | cout<<endl; | |
276 | ||
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; | |
281 | ||
0717295b | 282 | ef->Close(); |
3f7a17bb | 283 | TFile fc("AliESDComparison.root","RECREATE"); |
284 | c1->Write(); | |
285 | fc.Close(); | |
0717295b | 286 | |
4fa1e629 | 287 | gBenchmark->Stop("AliESDComparison"); |
288 | gBenchmark->Show("AliESDComparison"); | |
0717295b | 289 | |
4fa1e629 | 290 | delete rl; |
291 | return 0; | |
0717295b | 292 | } |