Macro to check the quality of ITS clusters
[u/mrichter/AliRoot.git] / ITS / ShowITSRecPoints.C
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
19 Int_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