.so cleanup: removed from gSystem->Load()
[u/mrichter/AliRoot.git] / ITS / UPGRADE / macros / efficiency.C
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
26 Int_t minNpoints=3;
27 Bool_t IsTrackable(TClonesArray *trackRef);
28 TArrayF radii;
29
30 void efficiency(char *title="",Bool_t isPrimary=kTRUE){
31
32   gSystem->Load("libITSUpgradeBase");
33   gSystem->Load("libITSUpgradeRec");
34   gSystem->Load("libITSUpgradeSim");
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;
174       if(track->GetNcls(0)<minNpoints) continue; // selection compatible with IsSelected
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 //___________________________________________________
377 Bool_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());
389     for(Int_t iLay=0;iLay<radii.GetSize();iLay++){
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
398   if(nTrackRef>=minNpoints) isOk=kTRUE;
399   return isOk; 
400 }
401