]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/UPGRADE/macros/efficiency.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Correcting Digits2Raw + prints commented
[u/mrichter/AliRoot.git] / ITS / UPGRADE / macros / efficiency.C
CommitLineData
13bc44df 1
2#if !defined(__CINT__) || defined(__MAKECINT__)
3#include <TH1F.h>
4#include <TH1I.h>
5#include <TH2D.h>
6#include <TFile.h>
7#include <Riostream.h>
8#include <TParticle.h>
9#include <TCanvas.h>
10#include <TLegend.h>
11#include <TSystem.h>
12#include <TROOT.h>
13#include <TClonesArray.h>
14#include <TTree.h>
15#include <TArrayF.h>
16#include <TStyle.h>
17#include "AliRun.h"
18#include "AliRunLoader.h"
19#include "AliStack.h"
20#include "AliESDEvent.h"
21#include "AliESDtrack.h"
22#include "AliTrackReference.h"
23#include "AliITSsegmentationUpgrade.h"
24#endif
25
26Bool_t IsTrackable(TClonesArray *trackRef);
27TArrayF radii;
28
29void efficiency(char *title="",Bool_t isPrimary=kTRUE){
30
31 gSystem->Load("libITSUpgradeBase.so");
32 gSystem->Load("libITSUpgradeRec.so");
33 gSystem->Load("libITSUpgradeSim.so");
34
35 // setting up the configuration to be considered
36 AliITSsegmentationUpgrade *seg = new AliITSsegmentationUpgrade();
37 if(!seg){
38 printf("no segmentation info available... Exiting");
39 return;
40 }
41 radii.Set(seg->GetNLayers());
42 for(Int_t l=0; l<seg->GetNLayers(); l++){
43 radii.AddAt(seg->GetRadius(l),l);
44 }
45 delete seg;
46
47 TString titlet(title);
48 titlet.ReplaceAll(" ","");
49
50 gROOT->SetStyle("Plain");
51 gStyle->SetPalette(1);
52 gStyle->SetOptStat(0);
53 gStyle->SetOptTitle(0);
54
55
56 //booking plots
57 // pt
58 Int_t nbins= 64;
59 Double_t ptmax = 1.6;
60
61 TH1F *hPtRec = new TH1F("hPtRec"," p_{T} distribution of reconstructed tracks",nbins,0,ptmax);
62 hPtRec->SetXTitle("p_{T}");
63 hPtRec->SetYTitle(Form("entries / %3.0f MeV/c",1000*(ptmax/nbins)));
64
65 TH1F *hPtRecMc = new TH1F("hPtRecMc"," MC p_{T} distribution of reconstructed tracks",nbins,0,ptmax);
66 hPtRecMc->SetXTitle("p_{T} GeV/c");
67 hPtRecMc->SetYTitle(Form("entries / %3.0f MeV/c",1000*(ptmax/nbins)));
68
69 TH1F *hPtRecGood = new TH1F("hPtRecGood"," P_{t} distribution of reconstructed tracks [positive labels]",nbins,0,ptmax);
70 hPtRecGood->SetXTitle("p_{T} GeV/c");
71 hPtRecGood->SetYTitle(Form("entries / %3.0f MeV/c",1000*(ptmax/nbins)));
72
73 TH1F *hPtMc = new TH1F("hPtMc", " p_{T} distribution of MC tracks",nbins,0,ptmax);
74 hPtMc->SetXTitle("p_{T} GeV/c");
75 hPtMc->SetYTitle(Form("entries / %3.0f MeV/c",1000*(ptmax/nbins)));
76
77 // Eta
78 Int_t nbinsEta = 24;
79 Double_t eta[2]={-1.2,1.2};
80
81 TH1F *hEtaRec = new TH1F("hEtaRec"," #eta distribution of reconstructed tracks",nbinsEta,eta[0],eta[1]);
82 hEtaRec->SetXTitle("#eta");
83 hEtaRec->SetYTitle(Form("entries / %3.3f ",((eta[1]-eta[0])/nbinsEta)));
84
85 TH1F *hEtaRecMc = new TH1F("hEtaRecMc"," #eta Mc distribution of reconstructed tracks",nbinsEta,eta[0],eta[1]);
86 hEtaRecMc->SetXTitle("#eta");
87 hEtaRecMc->SetYTitle(Form("entries / %3.3f ",((eta[1]-eta[0])/nbinsEta)));
88
89 TH1F *hEtaRecGood = new TH1F("hEtaRecGood"," #eta distribution of Good reconstructed tracks",nbinsEta,eta[0],eta[1]);
90 hEtaRecGood->SetXTitle("#eta");
91 hEtaRecGood->SetYTitle(Form("entries / %3.3f ",((eta[1]-eta[0])/nbinsEta)));
92
93 TH1F *hEtaMc = new TH1F("hEtaMc", " #eta distribution of MC tracks",nbinsEta,eta[0],eta[1]);
94 hEtaMc->SetXTitle("#eta");
95 hEtaMc->SetYTitle(Form("entries / %3.3f ",((eta[1]-eta[0])/nbinsEta)));
96
97
98 Double_t ptlimit=0.0;
99
100 Double_t cont=0;
101 Double_t ntrack=0;
102
103 // Init simulation
104 gAlice=NULL;
105 AliRunLoader* runLoader = AliRunLoader::Open("galice.root");
106 if (!runLoader) {
107 Error("OpengAlice", "no files, please check the path");
108 return ;
109 }
110 runLoader->LoadgAlice();
111 gAlice = runLoader->GetAliRun();
112 runLoader->LoadHeader();
113 runLoader->LoadKinematics();
114 runLoader->LoadTrackRefs();
115
116 //Trackref
117 TTree *trackRefTree = 0x0;
118 TClonesArray *trackRef = new TClonesArray("AliTrackReference",1000);
119
120 // ESD
121 TFile* esdFile = TFile::Open("AliESDs.root");
122 if (!esdFile || !esdFile->IsOpen()) {
123 Error("CheckESD", "opening ESD file %s failed", "AliESDs.root");
124 return ;
125 }
126 AliESDEvent * esd = new AliESDEvent;
127 TTree* tree = (TTree*) esdFile->Get("esdTree");
128 if (!tree) {
129 Error("CheckESD", "no ESD tree found");
130 return ;
131 }
132 esd->ReadFromTree(tree);
133
134 printf("setting the low pt value to %f \n",ptlimit);
135 // Loop over events
136
137 for(Int_t i=0; i<runLoader->GetNumberOfEvents(); i++){
138 //for(Int_t i=0; i< 1; i++){
139 cout << "Event "<< i << endl;
140 runLoader->GetEvent(i);
141
142 AliStack *stack = runLoader->Stack();
143 trackRefTree=runLoader->TreeTR();
144 TBranch *br = trackRefTree->GetBranch("TrackReferences");
145 if(!br) {
146 printf("no TR branch available , exiting \n");
147 return;
148 }
149 br->SetAddress(&trackRef);
150
151 Int_t nParticles = stack->GetNtrack();
152 trackRefTree=runLoader->TreeTR();
153 trackRefTree->SetBranchAddress("TrackReferences",&trackRef);
154 for(Int_t p=0; p<nParticles; p++){
155 if(stack->Particle(p)->Pt()<ptlimit) continue;
156 if(isPrimary && !stack->IsPhysicalPrimary(p)) continue;
157 trackRefTree->GetEntry(stack->TreeKEntry(p));
158
159 if(IsTrackable(trackRef)) {
160 cont++;
161 hPtMc->Fill((stack->Particle(p))->Pt());
162 hEtaMc->Fill((stack->Particle(p))->Eta());
163 }
164 }//loop on kine particle
165 cout << " # trackable "<< cont;
166
167 // ESD
168 tree->GetEvent(i);
169 for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
170 AliESDtrack* track = esd->GetTrack(iTrack);
171 if(track->Pt()<ptlimit) continue;
172 if(isPrimary && !stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) continue;
173 ntrack++;
174 hPtRec->Fill(track->Pt());
175 hEtaRec->Fill(track->Eta());
176 if(track->GetLabel()>=0) {
177 hPtRecGood->Fill(stack->Particle(track->GetLabel())->Pt());
178 hEtaRecGood->Fill(stack->Particle(track->GetLabel())->Eta());
179 }
180 hPtRecMc->Fill(stack->Particle(TMath::Abs(track->GetLabel()))->Pt());
181 hEtaRecMc->Fill(stack->Particle(TMath::Abs(track->GetLabel()))->Eta());
182 }// loop on tracks
183 cout << " tracks " << ntrack << endl;
184
185 }//loop events
186
187 cout << "Global Efficiency "<< ntrack/cont << endl;
188
189 TCanvas *c = new TCanvas();
190 c->cd();
191 c->cd()->SetLogy();
192 c->cd()->SetGridy();
193 c->cd()->SetGridx();
194 c->cd()->SetObjectStat(0);
195 hPtMc->SetMarkerStyle(20);
196 hPtMc->Draw("err");
197 hPtRec->SetLineColor(kRed);
198 hPtRec->SetMarkerStyle(20);
199 hPtRec->SetMarkerColor(kRed);
200 hPtRec->Sumw2();
201 hPtRec->Draw("sameerr");
202 hPtRecMc->SetLineColor(kBlue);
203 hPtRecMc->SetMarkerStyle(20);
204 hPtRecMc->SetMarkerColor(kBlue);
205 hPtRecMc->Sumw2();
206 hPtRecMc->Draw("sameerr");
207 hPtRecGood->SetLineColor(kGreen);
208 hPtRecGood->SetMarkerStyle(20);
209 hPtRecGood->SetMarkerColor(kGreen);
210 hPtRecGood->Sumw2();
211 hPtRecGood->Draw("sameerr");
212
213 TLegend *leg = new TLegend(0.5,0.7,0.99,0.95);
214 leg->SetHeader(Form("%s",title));
215 leg->AddEntry(hPtRecGood,"Good tracks (positive labels)","p");
216 leg->AddEntry(hPtRec,"Reconstructed tracks (Rec p_{T})","p");
217 leg->AddEntry(hPtRecMc,"Reconstructed tracks ","p");
218 leg->AddEntry(hPtMc,"MC trackable particles","p");
219 leg->Draw();
220 c->SaveAs(Form("ptdist_%s.png",titlet.Data()));
221
222 TCanvas *cc = new TCanvas();
223 cc->cd();
224 cc->cd()->SetGridy();
225 cc->cd()->SetGridx();
226
227 Double_t max = 1.2;
228
229 TH1F *effrec=new TH1F("hEffRec","tracking efficiency of reconstructed tracks",nbins,0,ptmax);
230 effrec->Divide(hPtRecMc,hPtMc,1,1,"b");
231 effrec->SetXTitle("p_{T} GeV/c");
232 effrec->SetLineColor(kBlue);
233 effrec->SetMarkerStyle(20);
234 effrec->SetMarkerColor(kBlue);
235 effrec->SetMaximum(max);
236 effrec->SetMinimum(0);
237 effrec->DrawCopy("err");
238
239
240 TH1F *effrecgood=new TH1F("hEffRecGood","tracking efficiency of reconstructed tracks [positive labels]",nbins,0,ptmax);
241 effrecgood->SetXTitle("P_{T}");
242 effrecgood->Divide(hPtRecGood,hPtMc,1,1,"b");
243 effrecgood->SetLineColor(kGreen);
244 effrecgood->SetMarkerStyle(20);
245 effrecgood->SetMarkerColor(kGreen);
246 effrecgood->SetMaximum(max);
247 effrecgood->SetMinimum(0);
248 effrecgood->DrawCopy("sameerr");
249
250 TH1F *purity=new TH1F("purity","track purity [positive labels / all tracks]",nbins,0,ptmax);
251 purity->Divide(hPtRecGood,hPtRecMc,1,1,"b");
252 purity->SetLineColor(kGray);
253 purity->SetMarkerStyle(28);
254 purity->SetMarkerColor(kOrange);
255 purity->SetMaximum(max);
256 purity->SetMinimum(0);
257 purity->DrawCopy("sameerr");
258
259 TH1F *fakes=new TH1F("fakes","fake ratio [negative labels / reconstructable]",nbins,0,ptmax);
260
261 fakes->Add(hPtRecMc,hPtRecGood,1,-1);
262 fakes->Divide(hPtMc);
263 fakes->SetLineColor(kViolet);
264 fakes->SetLineStyle(3);
265 fakes->SetLineWidth(3);
266 fakes->SetMarkerColor(kViolet);
267 fakes->SetMaximum(max);
268 fakes->SetMinimum(0);
269 fakes->DrawCopy("samehist");
270
271 TLegend *legEff = new TLegend(0.1,0.83,0.65,0.99);
272 legEff->SetHeader(Form("%s",title));
273 legEff->AddEntry(effrec,"Overall Efficiency ","p");
274 legEff->AddEntry(effrecgood,"Efficiency good tracks (positive labels)","p");
275 legEff->AddEntry(purity,"Purity (positive labels / All tracks)","p");
276 legEff->AddEntry(fakes,"Fakes","l");
277 legEff->Draw();
278
279 cc->SaveAs(Form("efficiency_%s.png",titlet.Data()));
280
281
282 TCanvas *cEta = new TCanvas();
283 cEta->cd();
284 cEta->cd()->SetGridy();
285 cEta->cd()->SetGridx();
286 cEta->cd()->SetObjectStat(0);
287 hEtaMc->SetMarkerStyle(20);
288 hEtaMc->Draw("err");
289 hEtaMc->SetMinimum(0);
290 hEtaRec->SetLineColor(kRed);
291 hEtaRec->SetMarkerStyle(20);
292 hEtaRec->SetMarkerColor(kRed);
293 hEtaRec->Sumw2();
294 hEtaRec->Draw("sameerr");
295 hEtaRecMc->SetLineColor(kBlue);
296 hEtaRecMc->SetMarkerStyle(20);
297 hEtaRecMc->SetMarkerColor(kBlue);
298 hEtaRecMc->Sumw2();
299 hEtaRecMc->Draw("sameerr");
300 hEtaRecGood->SetLineColor(kGreen);
301 hEtaRecGood->SetMarkerStyle(20);
302 hEtaRecGood->SetMarkerColor(kGreen);
303 hEtaRecGood->Sumw2();
304 hEtaRecGood->Draw("sameerr");
305
306 TLegend *legEta = new TLegend(0.2,0.3,0.59,0.55);
307 legEta->SetHeader(Form("%s",title));
308 legEta->AddEntry(hEtaRecGood,"Good tracks (positive labels)","p");
309 legEta->AddEntry(hEtaRec,"Reconstructed tracks (Rec p_{T})","p");
310 legEta->AddEntry(hEtaRecMc,"Reconstructed tracks ","p");
311 legEta->AddEntry(hEtaMc,"MC trackable particles","p");
312 legEta->Draw();
313 cEta->SaveAs(Form("etadist_%s.png",title));
314
315 TCanvas *ccEta = new TCanvas();
316 ccEta->cd();
317 ccEta->cd()->SetGridy();
318 ccEta->cd()->SetGridx();
319
320 TH1F *effrecEta=new TH1F("hEffRecEta","#eta tracking efficiency of reconstructed tracks",nbinsEta,eta[0],eta[1]);
321 effrecEta->Divide(hEtaRecMc,hEtaMc,1,1,"b");
322 effrecEta->SetXTitle("#eta");
323 effrecEta->SetLineColor(kBlue);
324 effrecEta->SetMarkerStyle(20);
325 effrecEta->SetMarkerColor(kBlue);
326 effrecEta->SetMaximum(max);
327 effrecEta->SetMinimum(0);
328 effrecEta->DrawCopy("err");
329
330
331 TH1F *effrecgoodEta=new TH1F("hEffRecGoodEta","#eta tracking efficiency of reconstructed tracks [positive labels]",nbinsEta,eta[0],eta[1]);
332 effrecgoodEta->SetXTitle("eta");
333 effrecgoodEta->Divide(hEtaRecGood,hEtaMc,1,1,"b");
334 effrecgoodEta->SetLineColor(kGreen);
335 effrecgoodEta->SetMarkerStyle(20);
336 effrecgoodEta->SetMarkerColor(kGreen);
337 effrecgoodEta->SetMaximum(max);
338 effrecgoodEta->SetMinimum(0);
339 effrecgoodEta->DrawCopy("sameerr");
340
341 TH1F *purityEta=new TH1F("purityEta","#eta track purity [positive labels / all tracks]",nbinsEta,eta[0],eta[1]);
342 purityEta->Divide(hEtaRecGood,hEtaRecMc,1,1,"b");
343 purityEta->SetLineColor(kGray);
344 purityEta->SetMarkerStyle(28);
345 purityEta->SetMarkerColor(kOrange);
346 purityEta->SetMaximum(max);
347 purityEta->SetMinimum(0);
348 purityEta->DrawCopy("sameerr");
349
350 TH1F *fakesEta=new TH1F("fakesEta","#eta fake ratio [negative labels / reconstructable]",nbinsEta,eta[0],eta[1]);
351 fakesEta->Add(hEtaRecMc,hEtaRecGood,1,-1);
352 fakesEta->Divide(hEtaMc);
353 fakesEta->SetLineColor(kViolet);
354 fakesEta->SetLineStyle(3);
355 fakesEta->SetLineWidth(3);
356 fakesEta->SetMarkerColor(kViolet);
357 fakesEta->SetMaximum(max);
358 fakesEta->SetMinimum(0);
359 fakesEta->DrawCopy("samehist");
360
361 TLegend *legEffEta = new TLegend(0.3,0.3,0.75,0.55);
362 legEffEta->SetHeader(Form("%s",title));
363 legEffEta->AddEntry(effrecEta,"Overall Efficiency ","p");
364 legEffEta->AddEntry(effrecgoodEta,"Efficiency good tracks (positive labels)","p");
365 legEffEta->AddEntry(purityEta,"Purity (positive labels / All tracks)","p");
366 legEffEta->AddEntry(fakesEta,"Fakes","l");
367 legEffEta->Draw();
368
369 ccEta->SaveAs(Form("efficiencyEta_%s.png",titlet.Data()));
370
371
372}// main
373
374//___________________________________________________
375Bool_t IsTrackable(TClonesArray *trackRef){
376
377 Bool_t isOk=kFALSE;
378
379 Int_t nTrackRef =0;
380 TArrayF isInLayer(radii.GetSize());
381 for(Int_t l=0; l<radii.GetSize(); l++ ) isInLayer.AddAt(0,l);
382 for(Int_t nT =0; nT<trackRef->GetEntries(); nT++){
383 AliTrackReference *trR = (AliTrackReference*)trackRef->At(nT);
384 if(!trR) continue;
385 if(trR->DetectorId()!=0)continue;
386 Double_t rPart = TMath::Sqrt(trR->X()*trR->X()+trR->Y()*trR->Y());
28e46359 387 for(Int_t iLay=0;iLay<radii.GetSize();iLay++){
13bc44df 388 if(TMath::Abs(rPart-radii.At(iLay))<0.01) isInLayer.AddAt(1,iLay);
389 }
390 }
391
392 for(Int_t iLayer=0;iLayer<radii.GetSize();iLayer++){
393 if(isInLayer.At(iLayer)) nTrackRef++;
394 }
395
396 if(nTrackRef>2) isOk=kTRUE;
397 return isOk;
398}
399