update for the NUA
[u/mrichter/AliRoot.git] / ITS / PlotITSHits.C
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
49 void 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 }