Converting PWG/TRD to native cmake
[u/mrichter/AliRoot.git] / ITS / PlotITSHits.C
CommitLineData
1603d531 1#if !defined(__CINT__) || defined(__MAKECINT__)
2#include<Riostream.h>
3#include<TROOT.h>
4#include<TArrayI.h>
5#include<TBranch.h>
6#include<TCanvas.h>
7#include<TClassTable.h>
8#include<TClonesArray.h>
9#include<TFile.h>
10#include<TStyle.h>
11#include<TH1.h>
12#include<TH2.h>
13#include<TLatex.h>
14#include <TInterpreter.h>
15#include <TGeoManager.h>
16#include<TObject.h>
17#include<TObjArray.h>
18#include<TTree.h>
19#include<TNtuple.h>
20#include<TParticle.h>
21#include "AliStack.h"
22#include "AliRun.h"
23#include "AliCDBManager.h"
24#include "AliGeomManager.h"
25#include "AliITS.h"
26#include "AliITSgeomTGeo.h"
27#include "AliITSDetTypeRec.h"
28#include "AliITSRecPoint.h"
29#include "AliITSRecPoint.h"
30#include "AliITSdigit.h"
31#include "AliITSdigitSSD.h"
32#include "AliITShit.h"
33#include "AliITSmodule.h"
34#include "AliITSsegmentation.h"
35#include "AliITSsegmentationSPD.h"
36#include "AliITSsegmentationSDD.h"
37#include "AliITSsegmentationSSD.h"
38#include "AliRunLoader.h"
39#include "AliITSLoader.h"
40#include "AliHeader.h"
41#endif
42
43/* $Id$ */
44
45// macro to display the coordindates (local+global) and the energy deposit
46// of ITS hits
47
48
49void PlotITSHits() {
50
51
52 AliCDBManager* man = AliCDBManager::Instance();
53 if (!man->IsDefaultStorageSet()) {
54 printf("Setting a local default storage and run number 0\n");
55 man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
56 man->SetRun(0);
57 }
58 else {
59 printf("Using deafult storage \n");
60 }
61 // retrives geometry
62 if(!gGeoManager){
63 AliGeomManager::LoadGeometry("geometry.root");
64 }
65 AliGeomManager::ApplyAlignObjsFromCDB("ITS");
66
67 Int_t totmod=AliITSgeomTGeo::GetNModules();
68 Int_t modmin=AliITSgeomTGeo::GetModuleIndex(1,1,1);
69 Int_t modmax=AliITSgeomTGeo::GetNModules()-1;
70
71 Float_t xlim[6]={4.5,7.5,16.,26.,40.,45.};
72 Float_t zlim[6]={15.,15.,22.,30.,45.,55.};
73
74 TH1F* hlayer=new TH1F("hlayer","",6,-0.5,5.5);
75 TH1F** hmod=new TH1F*[6];
76 TH1F** hxl=new TH1F*[6];
77 TH1F** hzl=new TH1F*[6];
78 TH1F** hxg=new TH1F*[6];
79 TH1F** hyg=new TH1F*[6];
80 TH1F** hzg=new TH1F*[6];
81 TH1F** hr=new TH1F*[6];
82 TH1F** hphi=new TH1F*[6];
83 TH2F** hzphi=new TH2F*[6];
84 TH1F** hq=new TH1F*[6];
85
86 Char_t name[10];
87 for(Int_t iLay=0;iLay<6;iLay++){
88 sprintf(name,"hmod%d",iLay+1);
89 hmod[iLay]=new TH1F(name,"",totmod,modmin-0.5,modmax+0.5);
90 hmod[iLay]->GetXaxis()->SetTitle("Module");
91 hmod[iLay]->GetXaxis()->CenterTitle();
92 sprintf(name,"hxloc%d",iLay+1);
93 hxl[iLay]=new TH1F(name,"",100,-4.,4.);
94 hxl[iLay]->GetXaxis()->SetTitle("Xloc (cm)");
95 hxl[iLay]->GetXaxis()->CenterTitle();
96 sprintf(name,"hzloc%d",iLay+1);
97 hzl[iLay]=new TH1F(name,"",100,-4.,4.);
98 hzl[iLay]->GetXaxis()->SetTitle("Zloc (cm)");
99 hzl[iLay]->GetXaxis()->CenterTitle();
100 sprintf(name,"hxgl%d",iLay+1);
101 hxg[iLay]=new TH1F(name,"",100,-xlim[iLay],xlim[iLay]);
102 hxg[iLay]->GetXaxis()->SetTitle("Xglob (cm)");
103 hxg[iLay]->GetXaxis()->CenterTitle();
104 sprintf(name,"hygl%d",iLay+1);
105 hyg[iLay]=new TH1F(name,"",100,-xlim[iLay],xlim[iLay]);
106 hyg[iLay]->GetXaxis()->SetTitle("Yglob (cm)");
107 hyg[iLay]->GetXaxis()->CenterTitle();
108 sprintf(name,"hzgl%d",iLay+1);
109 hzg[iLay]=new TH1F(name,"",100,-zlim[iLay],zlim[iLay]);
110 hzg[iLay]->GetXaxis()->SetTitle("Zglob (cm)");
111 hzg[iLay]->GetXaxis()->CenterTitle();
112 sprintf(name,"hr%d",iLay+1);
113 hr[iLay]=new TH1F(name,"",100,0.,50.);
114 hr[iLay]->GetXaxis()->SetTitle("r (cm)");
115 hr[iLay]->GetXaxis()->CenterTitle();
116 sprintf(name,"hphi%d",iLay+1);
117 hphi[iLay]=new TH1F(name,"",100,-TMath::Pi(),TMath::Pi());
118 hphi[iLay]->GetXaxis()->SetTitle("#varphi (rad)");
119 hphi[iLay]->GetXaxis()->CenterTitle();
120 sprintf(name,"hq%d",iLay+1);
121 hq[iLay]=new TH1F(name,"",100,0.,300.);
122 hq[iLay]->GetXaxis()->SetTitle("Charge (keV)");
123 hq[iLay]->GetXaxis()->CenterTitle();
124 sprintf(name,"hzphi%d",iLay+1);
125 hzphi[iLay]=new TH2F(name,Form("Layer %d",iLay+1),50,-TMath::Pi(),TMath::Pi(),50,-zlim[iLay],zlim[iLay]);
126 hzphi[iLay]->GetXaxis()->SetTitle("#varphi (rad)");
127 hzphi[iLay]->GetYaxis()->SetTitle("Zglob (cm)");
128 hzphi[iLay]->SetStats(0);
129 }
130
131
132 AliRunLoader* rl = AliRunLoader::Open("galice.root");
133 if (rl == 0x0){
134 cerr<<"Can not open session RL=NULL"<< endl;
135 return;
136 }
137 Int_t retval = rl->LoadgAlice();
138 if (retval){
139 cerr<<"LoadgAlice returned error"<<endl;
140 return;
141 }
142 gAlice=rl->GetAliRun();
143
144 retval = rl->LoadHeader();
145 if (retval){
146 cerr<<"LoadHeader returned error"<<endl;
147 return;
148 }
149
150 retval = rl->LoadKinematics();
151 if (retval){
152 cerr<<"LoadKinematics returned error"<<endl;
153 return;
154 }
155
156 AliITSLoader* ITSloader = (AliITSLoader*) rl->GetLoader("ITSLoader");
157 if(!ITSloader){
158 cerr<<"ITS loader not found"<<endl;
159 return;
160 }
161
162 ITSloader->LoadHits("read");
163 AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
164 ITS->SetTreeAddress();
165 Int_t totev=rl->GetNumberOfEvents();
166
167 for(Int_t iev=0; iev<totev; iev++){
168 rl->GetEvent(iev);
169
170 // HITS
171 TTree *TH = ITSloader->TreeH();
172 printf("Event %d Tracks %d\n",iev,(Int_t)TH->GetEntries());
173
174 // ITS
175 Int_t nmodules;
176 ITS->InitModules(-1,nmodules);
177 ITS->FillModules(0,0,nmodules," "," ");
178
179
180 for (Int_t mod=modmin; mod<=modmax; mod++){
181 Int_t lay,lad,det;
182 AliITSgeomTGeo::GetModuleId(mod,lay,lad,det);
183 lay--;
184
185 // Hits
186 AliITSmodule *modu = ITS->GetModule(mod);
187 TObjArray *arrHits = modu->GetHits();
188 Int_t nhits = arrHits->GetEntriesFast();
189 for (Int_t iHit=0;iHit<nhits;iHit++) {
190 AliITShit *hit = (AliITShit*) arrHits->At(iHit);
191 Int_t iMod=hit->GetModule();
192 hlayer->Fill(lay);
193 Double_t xl,yl,zl,xl0,yl0,zl0;
194 Double_t xg,yg,zg,tof,xg0,yg0,zg0,tof0;
195 Float_t hitloc[3],hitglo[3];
196 hit->GetPositionL(xl,yl,zl,tof);
197 hit->GetPositionL0(xl0,yl0,zl0,tof0);
198 hit->GetPositionG(xg,yg,zg,tof);
199 hit->GetPositionG0(xg0,yg0,zg0,tof0);
200 Double_t hitlen=TMath::Abs(yl-yl0);
201 if(hitlen<0.005) continue; // remove hits "shorter" than 50 um
202 if(lay > 1 && hitlen<0.025) continue; // remove hits "shorter" than 250 um in SDD,SSD
203 hitloc[0]=0.5*(xl+xl0);
204 hitloc[1]=0.5*(yl+yl0);
205 hitloc[2]=0.5*(zl+zl0);
206 hitglo[0]=0.5*(xg+xg0);
207 hitglo[1]=0.5*(yg+yg0);
208 hitglo[2]=0.5*(zg+zg0);
209 Float_t edep=hit->GetIonization()*1000000;
210 Float_t rad=TMath::Sqrt(hitglo[0]*hitglo[0]+hitglo[1]*hitglo[1]);
211 Float_t phi=TMath::ATan2(hitglo[1],hitglo[0]);
212 hmod[lay]->Fill(iMod);
213 hzl[lay]->Fill(hitloc[2]);
214 hxl[lay]->Fill(hitloc[0]);
215 hzg[lay]->Fill(hitglo[2]);
216 hyg[lay]->Fill(hitglo[1]);
217 hxg[lay]->Fill(hitglo[0]);
218 hr[lay]->Fill(rad);
219 hphi[lay]->Fill(phi);
220 hq[lay]->Fill(edep);
221 hzphi[lay]->Fill(phi,hitglo[2]);
222 }
223 }
224 }
225
226 gStyle->SetOptStat(10);
227 gStyle->SetPadBottomMargin(0.14);
228
229
230 TCanvas **c=new TCanvas*[6];
231 Char_t ctit[30];
232 for(Int_t iLay=0;iLay<6;iLay++){
233 sprintf(name,"can%d",iLay+1);
234 sprintf(ctit,"Layer %d",iLay+1);
235 c[iLay]=new TCanvas(name,ctit,1200,900);
236 c[iLay]->Divide(3,3,0.001,0.001);
237 c[iLay]->cd(1);
238 hmod[iLay]->Draw();
239 c[iLay]->cd(2);
240 hxl[iLay]->Draw();
241 c[iLay]->cd(3);
242 hzl[iLay]->Draw();
243 c[iLay]->cd(4);
244 hxg[iLay]->Draw();
245 c[iLay]->cd(5);
246 hyg[iLay]->Draw();
247 c[iLay]->cd(6);
248 hzg[iLay]->Draw();
249 c[iLay]->cd(7);
250 hr[iLay]->Draw();
251 c[iLay]->cd(8);
252 hphi[iLay]->Draw();
253 c[iLay]->cd(9);
254 hq[iLay]->Draw();
255 }
256
257 gStyle->SetPalette(1);
258 TLatex* tstat=new TLatex();
259 tstat->SetNDC();
260 TCanvas* cspd=new TCanvas("cspd","SPD",1000,600);
261 cspd->Divide(2,1);
262 cspd->cd(1);
263 hzphi[0]->Draw("colz");
264 tstat->DrawLatex(0.6,0.92,Form("# Clusters = %d",int(hzphi[0]->GetEntries())));
265 cspd->cd(2);
266 hzphi[1]->Draw("colz");
267 tstat->DrawLatex(0.6,0.92,Form("# Clusters = %d",int(hzphi[1]->GetEntries())));
268
269 TCanvas* csdd=new TCanvas("csdd","SDD",1000,600);
270 csdd->Divide(2,1);
271 csdd->cd(1);
272 hzphi[2]->Draw("colz");
273 tstat->DrawLatex(0.6,0.92,Form("# Clusters = %d",int(hzphi[2]->GetEntries())));
274 csdd->cd(2);
275 hzphi[3]->Draw("colz");
276 tstat->DrawLatex(0.6,0.92,Form("# Clusters = %d",int(hzphi[3]->GetEntries())));
277
278 TCanvas* cssd=new TCanvas("cssd","SSD",1000,600);
279 cssd->Divide(2,1);
280 cssd->cd(1);
281 hzphi[4]->Draw("colz");
282 tstat->DrawLatex(0.6,0.92,Form("# Clusters = %d",int(hzphi[4]->GetEntries())));
283 cssd->cd(2);
284 hzphi[5]->Draw("colz");
285 tstat->DrawLatex(0.6,0.92,Form("# Clusters = %d",int(hzphi[5]->GetEntries())));
286
287}