Updated directory names on LDC in SDD calibration plotting macro (Mario Sitta)
[u/mrichter/AliRoot.git] / ITS / macrosSDD / CheckSDDInESD.C
1 #if !defined(__CINT__) || defined(__MAKECINT__)
2 #include <TH1F.h>
3 #include <TH2F.h>
4 #include <TF1.h>
5 #include <TPad.h>
6 #include <TGraphErrors.h>
7 #include <TROOT.h>
8 #include <TFile.h>
9 #include <TTree.h>
10 #include <TMath.h>
11 #include <TCanvas.h>
12 #include <TStyle.h>
13 #include <TLatex.h>
14 #include "AliITSgeomTGeo.h"
15 #include "AliESDEvent.h"
16 #endif
17
18 enum {kAll, kTPCITS, kITSsa, kITSpureSA};
19
20 void CheckSDDInESD(TString filename="AliESDs.root", Int_t optTracks=kAll){
21
22
23   TFile* esdFile = TFile::Open(filename.Data());
24   if (!esdFile || !esdFile->IsOpen()) {
25     printf("Error in opening ESD file");
26     return;
27   }
28
29   AliESDEvent * esd = new AliESDEvent;
30   TTree* tree = (TTree*) esdFile->Get("esdTree");
31   if (!tree) {
32     printf("Error: no ESD tree found");
33     return;
34   }
35   esd->ReadFromTree(tree);
36   TH1F* hpt=new TH1F("hpt","",100,0.,10.);
37   TH1F* hphi=new TH1F("hphi","",100,-1,1);
38   TH1F* hlam=new TH1F("hlam","",100,-2.,2.);
39   TH1F* halpha=new TH1F("halpha","",100,-7,7);
40   TH1F* hitscl=new TH1F("hitscl","",7,-0.5,6.5);
41   TH1F* htpccl=new TH1F("htpccl","",200,-0.5,199.5);
42   TH1F* hitsmap=new TH1F("hitsmap","",64,-0.5,63.5);
43   TH1F* hclulay=new TH1F("hclulay","",7,-1.5,5.5);
44
45   TH1F* hvx=new TH1F("hvx","",100,-1.,1.);
46   TH1F* hvy=new TH1F("hvy","",100,-1.,1.);
47   TH1F* hvz=new TH1F("hvz","",100,-20.,20.);
48   TH1F* hdedx3=new TH1F("hdedx3","",100,0.,300.);
49   TH1F* hdedx4=new TH1F("hdedx4","",100,0.,300.);
50   TH1F* hdedx5=new TH1F("hdedx5","",100,0.,300.);
51   TH1F* hdedx6=new TH1F("hdedx6","",100,0.,300.);
52   TH1F* hStatus=new TH1F("hStatus","",11,-1.5,9.5);
53
54
55   // -- Local coordinates
56
57
58   // -- Module histos
59
60   TH1F* hAllPMod  = new TH1F("hAllPmod","Crossing Tracks vs. Module",260,239.5,499.5);
61   TH1F* hGoodPMod  = new TH1F("hGoodPmod","PointsAssocToTrack per Module",260,239.5,499.5);
62   TH1F* hBadRegMod  = new TH1F("hBadRegmod","Tracks in BadRegion per Module",260,239.5,499.5);
63   TH1F* hMissPMod  = new TH1F("hMissPmod","Missing Points per Module",260,239.5,499.5);
64   TH1F* hSkippedMod  = new TH1F("hSkippedmod","Tracks in Skipped Module",260,239.5,499.5);
65   TH1F* hOutAccMod  = new TH1F("hOutAccmod","Tracks outside zAcc per Module",260,239.5,499.5);
66   TH1F* hNoRefitMod  = new TH1F("hNoRefitmod","Points rejected in refit per Module",260,239.5,499.5);
67
68   TH1F* hAllPXloc  = new TH1F("hAllPxloc","Crossing Tracks vs. Xloc",75, -3.75, 3.75);
69   TH1F* hGoodPXloc  = new TH1F("hGoodPxloc","PointsAssocToTrack vs. Xloc",75, -3.75, 3.75);
70   TH1F* hBadRegXloc  = new TH1F("hBadRegxloc","Tracks in BadRegion vs. Xloc",75, -3.75, 3.75);
71   TH1F* hMissPXloc  = new TH1F("hMissPxloc","Missing Points vs. Xloc",75, -3.75, 3.75);
72   TH1F* hAllPZloc  = new TH1F("hAllPzloc","Crossing Tracks vs. Zloc",77, -3.85, 3.85);
73   TH1F* hGoodPZloc  = new TH1F("hGoodPzloc","PointsAssocToTrack vs. Zloc",77, -3.85, 3.85);
74   TH1F* hBadRegZloc  = new TH1F("hBadRegzloc","Tracks in BadRegion vs. Zloc",77, -3.85, 3.85);
75   TH1F* hMissPZloc  = new TH1F("hMissPzloc","Missing Points vs. Zloc",77, -3.85, 3.85);
76   TH2F* hdEdxVsMod=new TH2F("hdEdxVsMod","dE/dx vs. mod",260,239.5,499.5,100,0.,500.);
77
78   gStyle->SetPalette(1);
79   
80
81   for (Int_t iEvent = 0; iEvent < tree->GetEntries(); iEvent++) {
82     tree->GetEvent(iEvent);
83     if (!esd) {
84       printf("Error: no ESD object found for event %d", iEvent);
85       return;
86     }
87     cout<<"-------- Event "<<iEvent<<endl;
88     printf(" Tracks # = %d\n",esd->GetNumberOfTracks());
89     const AliESDVertex *spdv=esd->GetVertex();
90     printf(" SPD Primary Vertex in %f %f %f with %d contributors\n",spdv->GetXv(),spdv->GetYv(),spdv->GetZv(),spdv->GetNContributors());
91     const AliESDVertex *trkv=esd->GetPrimaryVertex();
92     printf(" Track Primary Vertex with %d contributors\n",trkv->GetNContributors());
93     if(spdv->IsFromVertexer3D()){
94       hvx->Fill(spdv->GetXv());
95       hvy->Fill(spdv->GetYv());
96       hvz->Fill(spdv->GetZv());
97     }
98     Double_t itss[4];
99     for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
100       AliESDtrack* track = esd->GetTrack(iTrack);
101       Int_t nITSclus=track->GetNcls(0);
102       UChar_t clumap=track->GetITSClusterMap();
103       Int_t nPointsForPid=0;
104       for(Int_t i=2; i<6; i++){
105         if(clumap&(1<<i)) ++nPointsForPid;
106       }
107       //      track->PropagateTo(4.,5.);
108       htpccl->Fill(track->GetNcls(1));
109       Int_t status=track->GetStatus();
110       Bool_t tpcin=0;
111       hStatus->Fill(-1.);
112       if(status & AliESDtrack::kTPCin){
113         tpcin=1;
114         hStatus->Fill(0.);
115       }
116       if(status & AliESDtrack::kTPCout){
117         hStatus->Fill(1.);
118       }
119       if(status & AliESDtrack::kTPCrefit){
120         hStatus->Fill(2.);
121       }
122       Bool_t itsin=0;
123       if(status & AliESDtrack::kITSin){
124         itsin=1;
125         hStatus->Fill(3.);
126       }
127       if(status & AliESDtrack::kITSout){
128         hStatus->Fill(4.);
129       }
130       if(status & AliESDtrack::kITSrefit){
131         hStatus->Fill(5.);
132       }
133       if(!tpcin && itsin){
134         hStatus->Fill(6.);
135       }
136       if(status & AliESDtrack::kITSpureSA){
137         hStatus->Fill(7.);
138       }
139
140       if(status & AliESDtrack::kITSrefit){
141         if((optTracks==kTPCITS) && !(status & AliESDtrack::kTPCin)) continue;
142         if((optTracks==kITSsa) && (status & AliESDtrack::kTPCin)) continue;
143         if((optTracks==kITSsa) && (status & AliESDtrack::kITSpureSA)) continue;
144         if((optTracks==kITSpureSA) && (status & AliESDtrack::kITSpureSA)) continue;
145
146          track->GetITSdEdxSamples(itss);
147         //      printf("Track %d (label %d) in ITS with %d clusters clumap %d pointspid= %d\n",iTrack,track->GetLabel(),nITSclus,clumap,nPointsForPid);
148         //printf("   dedx=%f %f %f %f\n",itss[0],itss[1],itss[2],itss[3]);
149         hitscl->Fill(nITSclus);
150         hdedx3->Fill(itss[0]);
151         hdedx4->Fill(itss[1]);
152         hdedx5->Fill(itss[2]);
153         hdedx6->Fill(itss[3]);
154         hitsmap->Fill(clumap);
155         hclulay->Fill(-1.);
156         for(Int_t iLay=0;iLay<6;iLay++){
157           if(clumap&1<<iLay) hclulay->Fill(iLay);
158         }
159         hpt->Fill(track->Pt());
160         hphi->Fill(TMath::ASin(track->GetSnp()));
161         hlam->Fill(TMath::ATan(track->GetTgl()));
162         halpha->Fill(track->GetAlpha());
163         Int_t iMod,status;
164         Float_t xloc,zloc;
165         for(Int_t iLay=2; iLay<=3; iLay++){
166           Bool_t ok=track->GetITSModuleIndexInfo(iLay,iMod,status,xloc,zloc);
167           if(ok){
168             iMod+=240;
169             hAllPMod->Fill(iMod);
170             hAllPXloc->Fill(xloc);
171             hAllPZloc->Fill(zloc);
172             if(status==1){
173               hGoodPMod->Fill(iMod);
174               hGoodPXloc->Fill(xloc);
175               hGoodPZloc->Fill(zloc);
176               if(track->Pt()>1.) hdEdxVsMod->Fill(iMod,itss[iLay-2]);
177             }
178             else if(status==2){ 
179               hBadRegMod->Fill(iMod);
180               hBadRegXloc->Fill(xloc);
181               hBadRegZloc->Fill(zloc);
182             }
183             else if(status==3) hSkippedMod->Fill(iMod);
184             else if(status==4) hOutAccMod->Fill(iMod);
185             else if(status==5){
186               hMissPMod->Fill(iMod);
187               hMissPXloc->Fill(xloc);
188               hMissPZloc->Fill(zloc);
189             }
190             else if(status==6) hNoRefitMod->Fill(iMod);
191           }
192         }
193       }
194     }
195   }
196   Float_t norm=hclulay->GetBinContent(1);
197   if(norm<1.) norm=1.;
198   hclulay->Scale(1./norm);
199   gStyle->SetLineWidth(2);
200
201   TCanvas* c1=new TCanvas("c1","Track quantities",900,900);
202   c1->Divide(2,2);
203   c1->cd(1);
204   htpccl->Draw();
205   htpccl->GetXaxis()->SetTitle("Clusters in TPC ");
206   c1->cd(2);
207   hitscl->Draw();
208   hitscl->GetXaxis()->SetTitle("Clusters in ITS ");
209   c1->cd(3);
210   hclulay->Draw();
211   hclulay->GetXaxis()->SetRange(2,7);
212   hclulay->GetXaxis()->SetTitle("# ITS Layer");
213   hclulay->GetYaxis()->SetTitle("Fraction of tracks with point in Layer x");
214   c1->cd(4);
215
216   TCanvas* c2=new TCanvas("c2","dedx per Layer",900,900);
217   c2->Divide(2,2);
218   c2->cd(1);
219   hdedx3->Draw();
220   hdedx3->GetXaxis()->SetTitle("dE/dx Lay3");
221   c2->cd(2);
222   hdedx4->Draw();
223   hdedx4->GetXaxis()->SetTitle("dE/dx Lay4");
224   c2->cd(3);
225   hdedx5->Draw();
226   hdedx5->GetXaxis()->SetTitle("dE/dx Lay5");
227   c2->cd(4);
228   hdedx6->Draw();
229   hdedx6->GetXaxis()->SetTitle("dE/dx Lay6");
230
231   hdEdxVsMod->SetStats(0);
232   TCanvas* cdedx=new TCanvas("cdedx","dedx SDD",1400,600);
233   cdedx->SetLogz();
234   hdEdxVsMod->Draw("col"); 
235   hdEdxVsMod->GetXaxis()->SetTitle("SDD Module Id");
236   hdEdxVsMod->GetYaxis()->SetTitle("dE/dx (keV/300 #mum)");
237   hdEdxVsMod->GetYaxis()->SetTitleOffset(1.25);
238
239
240
241   TCanvas* cv=new TCanvas("cv","Vertex",600,900);
242   cv->Divide(1,3);
243   cv->cd(1);
244   hvx->Draw();
245   hvx->GetXaxis()->SetTitle("Xv (cm)");
246   cv->cd(2);
247   hvy->Draw();
248   hvy->GetXaxis()->SetTitle("Yv (cm)");
249   cv->cd(3);
250   hvz->Draw();
251   hvz->GetXaxis()->SetTitle("Xv (cm)");
252
253   hGoodPMod->SetStats(0);
254   hGoodPMod->SetTitle("");
255   TCanvas* ceff0=new TCanvas("ceff0","ModuleIndexInfo",1000,600);
256   hGoodPMod->Draw("e");
257   hGoodPMod->GetXaxis()->SetTitle("SDD Module Id");
258   hGoodPMod->GetYaxis()->SetTitle("Number of tracks");
259   hMissPMod->SetLineColor(2);
260   hMissPMod->SetMarkerColor(2);
261   hMissPMod->SetMarkerStyle(22);
262   hMissPMod->SetMarkerSize(0.5);
263   hMissPMod->Draw("psame");
264   hBadRegMod->SetLineColor(kGreen+1);
265   hBadRegMod->SetMarkerColor(kGreen+1);
266   hBadRegMod->SetMarkerStyle(20);
267   hBadRegMod->SetMarkerSize(0.5);
268   hBadRegMod->Draw("esame");
269   hSkippedMod->SetLineColor(kYellow);
270   hSkippedMod->Draw("esame");
271   hOutAccMod->SetLineColor(4);
272   hOutAccMod->Draw("esame");
273   hNoRefitMod->SetLineColor(6);
274   hNoRefitMod->Draw("esame");
275   TLatex* t1=new TLatex(0.7,0.85,"Good Point");
276   t1->SetNDC();
277   t1->SetTextColor(1);
278   t1->Draw();
279   TLatex* t2=new TLatex(0.7,0.8,"Missing Point");
280   t2->SetNDC();
281   t2->SetTextColor(2);
282   t2->Draw();
283   TLatex* t3=new TLatex(0.7,0.75,"Bad Region");
284   t3->SetNDC();
285   t3->SetTextColor(kGreen+1);
286   t3->Draw();
287   ceff0->Update();
288
289   TH1F* heff=new TH1F("heff","",260,239.5,499.5);
290   for(Int_t imod=0; imod<260;imod++){
291     Float_t numer=hGoodPMod->GetBinContent(imod+1)+hBadRegMod->GetBinContent(imod+1)+hOutAccMod->GetBinContent(imod+1)+hNoRefitMod->GetBinContent(imod+1);
292     Float_t denom=hAllPMod->GetBinContent(imod+1);
293     Float_t eff=0.;
294     Float_t erreff=0.;
295     if(denom>0){
296       eff=numer/denom;
297       erreff=TMath::Sqrt(eff*(1-eff)/denom);
298     }
299     heff->SetBinContent(imod+1,eff);
300     heff->SetBinError(imod+1,erreff);
301   }
302
303   printf("---- Modules with efficiency < 90%% ----\n");
304   heff->SetStats(0);
305   TCanvas* ceff1=new TCanvas("ceff1","Efficiency",1000,600);
306   heff->Draw();
307   heff->GetXaxis()->SetTitle("SDD Module Id");
308   heff->GetYaxis()->SetTitle("Fraction of tracks with point in good region");
309   for(Int_t ibin=1; ibin<=heff->GetNbinsX(); ibin++){
310     Float_t e=heff->GetBinContent(ibin);
311     if(e<0.9){
312       Int_t iMod=(Int_t)heff->GetBinCenter(ibin);
313       Int_t lay,lad,det;
314       AliITSgeomTGeo::GetModuleId(iMod,lay,lad,det);
315       printf("Module %d - Layer %d Ladder %2d Det %d  -   Eff. %.3f\n",iMod,lay,lad,det,heff->GetBinContent(ibin));
316     }
317   }
318   ceff1->Update();
319
320
321
322   hGoodPXloc->SetTitle("");
323   hGoodPZloc->SetTitle("");
324   hGoodPXloc->SetStats(0);
325   hGoodPZloc->SetStats(0);
326   hGoodPXloc->SetMinimum(0);
327   hGoodPZloc->SetMinimum(0);
328   TCanvas* ceff2=new TCanvas("ceff2","LocalCoord",1000,600);
329   ceff2->Divide(2,1);
330   ceff2->cd(1);
331   hGoodPXloc->Draw("e");
332   hGoodPXloc->GetXaxis()->SetTitle("Xlocal (cm)");
333   hGoodPXloc->GetYaxis()->SetTitle("Number of tracks");
334   hMissPXloc->SetLineColor(2);
335   hMissPXloc->SetMarkerColor(2);
336   hMissPXloc->SetMarkerStyle(22);
337   hMissPXloc->SetMarkerSize(0.5);
338   hMissPXloc->Draw("psame");
339   hBadRegXloc->SetLineColor(kGreen+1);
340   hBadRegXloc->SetMarkerColor(kGreen+1);
341   hBadRegXloc->SetMarkerStyle(20);
342   hBadRegXloc->SetMarkerSize(0.5);
343   hBadRegXloc->Draw("psame");
344   t1->Draw();
345   t2->Draw();
346   t3->Draw();
347   ceff2->cd(2);
348   hGoodPZloc->Draw("e");
349   hGoodPZloc->GetXaxis()->SetTitle("Zlocal (cm)");
350   hGoodPZloc->GetYaxis()->SetTitle("Number of tracks");
351   hMissPZloc->SetLineColor(2);
352   hMissPZloc->SetMarkerColor(2);
353   hMissPZloc->SetMarkerStyle(22);
354   hMissPZloc->SetMarkerSize(0.5);
355   hMissPZloc->Draw("psame");
356   hBadRegZloc->SetLineColor(kGreen+1);
357   hBadRegZloc->SetMarkerColor(kGreen+1);
358   hBadRegZloc->SetMarkerStyle(20);
359   hBadRegZloc->SetMarkerSize(0.5);
360   hBadRegZloc->Draw("psame");
361   t1->Draw();
362   t2->Draw();
363   t3->Draw();
364   ceff2->Update();
365
366
367   
368 }
369
370
371
372
373