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