Macro to check the quality of ITS clusters
[u/mrichter/AliRoot.git] / ITS / ShowITSRecPoints.C
CommitLineData
29d63bb9 1#if !defined(__CINT__) || defined(__MAKECINT__)
2#include <TCanvas.h>
3#include <TClassTable.h>
4#include <TGraph.h>
5#include <TGraph2D.h>
6#include <TGeoManager.h>
7#include <TH1.h>
8#include <TInterpreter.h>
9#include "AliCDBManager.h"
10#include "AliGeomManager.h"
11#include "AliHeader.h"
12#include "AliITS.h"
13#include "AliITSDetTypeRec.h"
14#include "AliITSgeom.h"
15#include "AliITSRecPoint.h"
16#include "AliRun.h"
17#endif
18
19Int_t ShowITSRecPoints(){
20 ///////////////////////////////////////////////////////////////////////
21 // Macro to check clusters in the 6 ITS layers //
22 // Provides: //
23 // 6 canvases with 9 plots each (1 canvas for each layer) //
24 // 3 canvases with cluster XY coordinates for the first 3 events //
25 ///////////////////////////////////////////////////////////////////////
26
27 if (gClassTable->GetID("AliRun") < 0) {
28 gInterpreter->ExecuteMacro("loadlibs.C");
29 }
30 else {
31 if(gAlice){
32 delete gAlice->GetRunLoader();
33 delete gAlice;
34 gAlice=0;
35 }
36 }
37 // Set OCDB if needed
38 AliCDBManager* man = AliCDBManager::Instance();
39 if (!man->IsDefaultStorageSet()) {
40 printf("Setting a local default storage and run number 0\n");
41 man->SetDefaultStorage("local://$ALICE_ROOT");
42 man->SetRun(0);
43 }
44 else {
45 printf("Using deafult storage \n");
46 }
47
48 // retrives geometry
49 if(!gGeoManager){
50 AliGeomManager::LoadGeometry("geometry.root");
51 }
52
53 AliRunLoader* rl = AliRunLoader::Open("galice.root");
54 if (rl == 0x0){
55 cerr<<"Can not open session RL=NULL"<< endl;
56 return -1;
57 }
58 Int_t retval = rl->LoadgAlice();
59 if (retval){
60 cerr<<"LoadgAlice returned error"<<endl;
61 return -1;
62 }
63 gAlice=rl->GetAliRun();
64
65 retval = rl->LoadHeader();
66 if (retval){
67 cerr<<"LoadHeader returned error"<<endl;
68 return -1;
69 }
70
71
72 AliITSLoader* ITSloader = (AliITSLoader*) rl->GetLoader("ITSLoader");
73 if(!ITSloader){
74 cerr<<"ITS loader not found"<<endl;
75 return -1;
76 }
77 ITSloader->LoadRecPoints("read");
78
79 Float_t cluglo[3]={0.,0.,0.};
80 AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
81 ITS->SetTreeAddress();
82 AliITSgeom *geom = ITS->GetITSgeom();
83 AliITSDetTypeRec* detTypeRec = new AliITSDetTypeRec();
84 detTypeRec->SetITSgeom(ITSloader->GetITSgeom());
85 detTypeRec->SetDefaults();
86
87 Int_t modmin=geom->GetStartDet(0);
88 Int_t modmax=geom->GetLastDet(2);
89 Int_t totmod=modmax-modmin;
90 Float_t xlim[6]={4.5,7.5,16.,26.,40.,45.};
91 Float_t zlim[6]={15.,15.,22.,30.,45.,55.};
92
93 TH1F* hlayer=new TH1F("hlayer","",6,0.5,6.5);
94 TH1F** hmod=new TH1F*[6];
95 TH1F** hxl=new TH1F*[6];
96 TH1F** hzl=new TH1F*[6];
97 TH1F** hxg=new TH1F*[6];
98 TH1F** hyg=new TH1F*[6];
99 TH1F** hzg=new TH1F*[6];
100 TH1F** hr=new TH1F*[6];
101 TH1F** hphi=new TH1F*[6];
102 TH1F** hq=new TH1F*[6];
103
104 Char_t name[10];
105 for(Int_t iLay=0;iLay<6;iLay++){
106 sprintf(name,"hmod%d",iLay+1);
107 hmod[iLay]=new TH1F(name,"",totmod,modmin-0.5,modmax+0.5);
108 hmod[iLay]->GetXaxis()->SetTitle("Module");
109 hmod[iLay]->GetXaxis()->CenterTitle();
110 sprintf(name,"hxloc%d",iLay+1);
111 hxl[iLay]=new TH1F(name,"",100,-4.,4.);
112 hxl[iLay]->GetXaxis()->SetTitle("Xloc (cm)");
113 hxl[iLay]->GetXaxis()->CenterTitle();
114 sprintf(name,"hzloc%d",iLay+1);
115 hzl[iLay]=new TH1F(name,"",100,-4.,4.);
116 hzl[iLay]->GetXaxis()->SetTitle("Zloc (cm)");
117 hzl[iLay]->GetXaxis()->CenterTitle();
118 sprintf(name,"hxgl%d",iLay+1);
119 hxg[iLay]=new TH1F(name,"",100,-xlim[iLay],xlim[iLay]);
120 hxg[iLay]->GetXaxis()->SetTitle("Xglob (cm)");
121 hxg[iLay]->GetXaxis()->CenterTitle();
122 sprintf(name,"hygl%d",iLay+1);
123 hyg[iLay]=new TH1F(name,"",100,-xlim[iLay],xlim[iLay]);
124 hyg[iLay]->GetXaxis()->SetTitle("Yglob (cm)");
125 hyg[iLay]->GetXaxis()->CenterTitle();
126 sprintf(name,"hzgl%d",iLay+1);
127 hzg[iLay]=new TH1F(name,"",100,-zlim[iLay],zlim[iLay]);
128 hzg[iLay]->GetXaxis()->SetTitle("Zglob (cm)");
129 hzg[iLay]->GetXaxis()->CenterTitle();
130 sprintf(name,"hr%d",iLay+1);
131 hr[iLay]=new TH1F(name,"",100,0.,50.);
132 hr[iLay]->GetXaxis()->SetTitle("r (cm)");
133 hr[iLay]->GetXaxis()->CenterTitle();
134 sprintf(name,"hphi%d",iLay+1);
135 hphi[iLay]=new TH1F(name,"",100,-TMath::Pi(),TMath::Pi());
136 hphi[iLay]->GetXaxis()->SetTitle("#varphi (rad)");
137 hphi[iLay]->GetXaxis()->CenterTitle();
138 sprintf(name,"hq%d",iLay+1);
139 hq[iLay]=new TH1F(name,"",100,0.,300.);
140 hq[iLay]->GetXaxis()->SetTitle("Charge (keV)");
141 hq[iLay]->GetXaxis()->CenterTitle();
142 }
143
144 TGraph **gpts=new TGraph*[3];
145 TGraph2D **gpts3d=new TGraph2D*[3];
146 for(Int_t i=0;i<3;i++){
147 gpts[i]=new TGraph(0);
148 gpts3d[i]=new TGraph2D(0);
149 }
150 Int_t totev=rl->GetNumberOfEvents();
151 printf("Total Number of events = %d\n",totev);
152
153 for(Int_t iev=0;iev<totev;iev++){
154 rl->GetEvent(iev);
155 TTree *TR = ITSloader->TreeR();
156 TClonesArray *ITSrec = detTypeRec->RecPoints();
157 TBranch *branch = 0;
158 if(TR && ITSrec){
159 branch = ITSloader->TreeR()->GetBranch("ITSRecPoints");
160 if(branch)branch->SetAddress(&ITSrec);
161 }
162 Int_t nparticles = rl->GetHeader()->GetNtrack();
163 cout<<"Event #"<<iev<<" #Particles="<<nparticles<<endl;
164
165
166 Int_t ipt=0;
167 for(Int_t subd=0;subd<3;subd++){
168
169 Int_t first = geom->GetStartDet(subd);
170 Int_t last = geom->GetLastDet(subd);
171
172 for (Int_t mod=first; mod<=last; mod++){
173 detTypeRec->ResetRecPoints();
174 branch->GetEvent(mod);
175 Int_t nrecp = ITSrec->GetEntries();
176 if(nrecp>0){
177 for(Int_t irec=0;irec<nrecp;irec++) {
178 AliITSRecPoint *recp = (AliITSRecPoint*)ITSrec->At(irec);
179 Int_t lay=recp->GetLayer();
180 hlayer->Fill(lay);
181 recp->GetGlobalXYZ(cluglo);
182 Float_t rad=TMath::Sqrt(cluglo[0]*cluglo[0]+cluglo[1]*cluglo[1]);
183 Float_t phi=TMath::ATan2(cluglo[1],cluglo[0]);
184 if(iev<3){
185 gpts[iev]->SetPoint(ipt,cluglo[0],cluglo[1]);
186 gpts3d[iev]->SetPoint(ipt,cluglo[0],cluglo[1],cluglo[2]);
187 ipt++;
188 }
189 hmod[lay]->Fill(mod);
190 hzl[lay]->Fill(recp->GetDetLocalZ());
191 hxl[lay]->Fill(recp->GetDetLocalX());
192 hzg[lay]->Fill(cluglo[2]);
193 hyg[lay]->Fill(cluglo[1]);
194 hxg[lay]->Fill(cluglo[0]);
195 hr[lay]->Fill(rad);
196 hphi[lay]->Fill(phi);
197 hq[lay]->Fill(recp->GetQ());
198 }
199 }
200 }
201 }
202 }
203 TCanvas **c=new TCanvas*[6];
204 Char_t ctit[30];
205 for(Int_t iLay=0;iLay<6;iLay++){
206 sprintf(name,"can%d",iLay+1);
207 sprintf(ctit,"Layer %d",iLay+1);
208 c[iLay]=new TCanvas(name,ctit,1200,900);
209 c[iLay]->Divide(3,3,0.001,0.001);
210 c[iLay]->cd(1);
211 hmod[iLay]->Draw();
212 c[iLay]->cd(2);
213 hxl[iLay]->Draw();
214 c[iLay]->cd(3);
215 hzl[iLay]->Draw();
216 c[iLay]->cd(4);
217 hxg[iLay]->Draw();
218 c[iLay]->cd(5);
219 hyg[iLay]->Draw();
220 c[iLay]->cd(6);
221 hzg[iLay]->Draw();
222 c[iLay]->cd(7);
223 hr[iLay]->Draw();
224 c[iLay]->cd(8);
225 hphi[iLay]->Draw();
226 c[iLay]->cd(9);
227 hq[iLay]->Draw();
228 }
229
230 TCanvas *cev0;
231 cev0=new TCanvas("cev0","Event 0",600,600);
232 gpts[0]->SetMarkerStyle(7);
233 gpts[0]->SetTitle(0);
234 gpts[0]->Draw("AP");
235
236 TCanvas *cev1;
237 cev1=new TCanvas("cev1","Event 1",600,600);
238 gpts[1]->SetMarkerStyle(7);
239 gpts[1]->SetTitle(0);
240 gpts[1]->Draw("AP");
241
242 TCanvas *cev2;
243 cev2=new TCanvas("cev2","Event 2",600,600);
244 gpts[2]->SetMarkerStyle(7);
245 gpts[2]->SetTitle(0);
246 gpts[2]->Draw("AP");
247
248
249 return 0;
250}
251
252