Replace AliITSgeom with AliITSgeomTGeo. Possibility to select the event to be displa...
[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 "AliGeomManager.h"
10 #include "AliHeader.h"
11 #include "AliITS.h"
12 #include "AliITSDetTypeRec.h"
13 #include "AliITSgeomTGeo.h"
14 #include "AliITSRecPoint.h"
15 #include "AliRun.h"
16 #endif
17
18 Int_t ShowITSRecPoints(Int_t nevfordisp=0){
19   ///////////////////////////////////////////////////////////////////////
20   // Macro to check clusters in the 6 ITS layers                       //
21   // Provides:                                                         //
22   //  6 canvases with 9 plots each (1 canvas for each layer)           //
23   //  3 canvases with cluster XY coordinates for the first 3 events    //
24   ///////////////////////////////////////////////////////////////////////
25
26   if (gClassTable->GetID("AliRun") < 0) {
27     gInterpreter->ExecuteMacro("loadlibs.C");
28   }
29  
30   // retrives geometry 
31   if(!gGeoManager){
32     AliGeomManager::LoadGeometry("geometry.root");
33   }
34
35   AliRunLoader* rl = AliRunLoader::Open("galice.root");
36   if (rl == 0x0){
37     cerr<<"Can not open session RL=NULL"<< endl;
38     return -1;
39   }
40
41   AliITSLoader* ITSloader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
42   if(!ITSloader){
43     cerr<<"ITS loader not found"<<endl;
44     return -1;
45   }
46   ITSloader->LoadRecPoints("read");
47
48   Float_t cluglo[3]={0.,0.,0.}; 
49   AliITSDetTypeRec* detTypeRec = new AliITSDetTypeRec();
50
51   Int_t totmod=AliITSgeomTGeo::GetNModules();
52   Int_t modmin=AliITSgeomTGeo::GetModuleIndex(1,1,1);
53   Int_t modmax=AliITSgeomTGeo::GetNModules()-1;
54
55   Float_t xlim[6]={4.5,7.5,16.,26.,40.,45.};
56   Float_t zlim[6]={15.,15.,22.,30.,45.,55.};
57
58   TH1F* hlayer=new TH1F("hlayer","",6,-0.5,5.5);
59   TH1F** hmod=new TH1F*[6];
60   TH1F** hxl=new TH1F*[6];
61   TH1F** hzl=new TH1F*[6];
62   TH1F** hxg=new TH1F*[6];
63   TH1F** hyg=new TH1F*[6];
64   TH1F** hzg=new TH1F*[6];
65   TH1F** hr=new TH1F*[6];
66   TH1F** hphi=new TH1F*[6];
67   TH1F** hq=new TH1F*[6];
68
69   Char_t name[10];
70   for(Int_t iLay=0;iLay<6;iLay++){
71     sprintf(name,"hmod%d",iLay+1);
72     hmod[iLay]=new TH1F(name,"",totmod,modmin-0.5,modmax+0.5);
73     hmod[iLay]->GetXaxis()->SetTitle("Module");
74     hmod[iLay]->GetXaxis()->CenterTitle();
75     sprintf(name,"hxloc%d",iLay+1);
76     hxl[iLay]=new TH1F(name,"",100,-4.,4.);
77     hxl[iLay]->GetXaxis()->SetTitle("Xloc (cm)");
78     hxl[iLay]->GetXaxis()->CenterTitle();
79     sprintf(name,"hzloc%d",iLay+1);
80     hzl[iLay]=new TH1F(name,"",100,-4.,4.);
81     hzl[iLay]->GetXaxis()->SetTitle("Zloc (cm)");
82     hzl[iLay]->GetXaxis()->CenterTitle();
83     sprintf(name,"hxgl%d",iLay+1);
84     hxg[iLay]=new TH1F(name,"",100,-xlim[iLay],xlim[iLay]);
85     hxg[iLay]->GetXaxis()->SetTitle("Xglob (cm)");
86     hxg[iLay]->GetXaxis()->CenterTitle();
87     sprintf(name,"hygl%d",iLay+1);
88     hyg[iLay]=new TH1F(name,"",100,-xlim[iLay],xlim[iLay]);
89     hyg[iLay]->GetXaxis()->SetTitle("Yglob (cm)");
90     hyg[iLay]->GetXaxis()->CenterTitle();
91     sprintf(name,"hzgl%d",iLay+1);
92     hzg[iLay]=new TH1F(name,"",100,-zlim[iLay],zlim[iLay]);
93     hzg[iLay]->GetXaxis()->SetTitle("Zglob (cm)");
94     hzg[iLay]->GetXaxis()->CenterTitle();
95     sprintf(name,"hr%d",iLay+1);
96     hr[iLay]=new TH1F(name,"",100,0.,50.);
97     hr[iLay]->GetXaxis()->SetTitle("r (cm)");
98     hr[iLay]->GetXaxis()->CenterTitle();
99     sprintf(name,"hphi%d",iLay+1);
100     hphi[iLay]=new TH1F(name,"",100,-TMath::Pi(),TMath::Pi());    
101     hphi[iLay]->GetXaxis()->SetTitle("#varphi (rad)");
102     hphi[iLay]->GetXaxis()->CenterTitle();
103     sprintf(name,"hq%d",iLay+1);
104     hq[iLay]=new TH1F(name,"",100,0.,300.);    
105     hq[iLay]->GetXaxis()->SetTitle("Charge (keV)");
106     hq[iLay]->GetXaxis()->CenterTitle();
107   }
108
109   TGraph *gptsXY=new TGraph(0);
110   TGraph *gptsRZ=new TGraph(0);
111
112   Int_t totev=rl->GetNumberOfEvents();
113   printf("Total Number of events = %d\n",totev);
114
115   for(Int_t iev=0;iev<totev;iev++){
116     rl->GetEvent(iev);
117     TTree *TR = ITSloader->TreeR();
118     TClonesArray *ITSrec  = detTypeRec->RecPoints();
119     TBranch *branch = 0;
120     if(TR && ITSrec){
121       branch = ITSloader->TreeR()->GetBranch("ITSRecPoints");
122       if(branch)branch->SetAddress(&ITSrec);
123     }
124     Int_t nparticles = rl->GetHeader()->GetNtrack();
125     cout<<"Event #"<<iev<<"   #Particles="<<nparticles<<endl;
126
127
128     Int_t ipt=0;
129     for (Int_t mod=modmin; mod<=modmax; mod++){
130       detTypeRec->ResetRecPoints();
131       branch->GetEvent(mod);
132       Int_t nrecp = ITSrec->GetEntries();
133       if(nrecp>0){
134         for(Int_t irec=0;irec<nrecp;irec++) {
135           AliITSRecPoint *recp = (AliITSRecPoint*)ITSrec->At(irec);
136           Int_t lay=recp->GetLayer();
137           hlayer->Fill(lay);
138           recp->GetGlobalXYZ(cluglo);
139           Float_t rad=TMath::Sqrt(cluglo[0]*cluglo[0]+cluglo[1]*cluglo[1]); 
140           Float_t phi=TMath::ATan2(cluglo[1],cluglo[0]);
141           if(iev==nevfordisp){
142             gptsXY->SetPoint(ipt,cluglo[0],cluglo[1]);
143             if(cluglo[1]>0) gptsRZ->SetPoint(ipt,cluglo[2],rad);
144             else gptsRZ->SetPoint(ipt,cluglo[2],-rad);
145             ipt++;
146           }
147           hmod[lay]->Fill(mod);
148           hzl[lay]->Fill(recp->GetDetLocalZ());
149           hxl[lay]->Fill(recp->GetDetLocalX());
150           hzg[lay]->Fill(cluglo[2]);
151           hyg[lay]->Fill(cluglo[1]);
152           hxg[lay]->Fill(cluglo[0]);
153           hr[lay]->Fill(rad);
154           hphi[lay]->Fill(phi);
155           hq[lay]->Fill(recp->GetQ());
156         }
157       }
158     }
159   }
160   
161   TCanvas **c=new TCanvas*[6];
162   Char_t ctit[30];
163   for(Int_t iLay=0;iLay<6;iLay++){
164     sprintf(name,"can%d",iLay+1);
165     sprintf(ctit,"Layer %d",iLay+1);
166     c[iLay]=new TCanvas(name,ctit,1200,900);
167     c[iLay]->Divide(3,3,0.001,0.001);
168     c[iLay]->cd(1);
169     hmod[iLay]->Draw();
170     c[iLay]->cd(2);
171     hxl[iLay]->Draw();
172     c[iLay]->cd(3);
173     hzl[iLay]->Draw();
174     c[iLay]->cd(4);
175     hxg[iLay]->Draw();
176     c[iLay]->cd(5);
177     hyg[iLay]->Draw();
178     c[iLay]->cd(6);
179     hzg[iLay]->Draw();
180     c[iLay]->cd(7);
181     hr[iLay]->Draw();
182     c[iLay]->cd(8);
183     hphi[iLay]->Draw();    
184     c[iLay]->cd(9);
185     hq[iLay]->Draw();    
186   }
187
188   TCanvas *cev0;
189   sprintf(ctit,"Event %d XY",nevfordisp);
190   cev0=new TCanvas("cev0",ctit,600,600);
191   if(gptsXY->GetN()>0){
192     gptsXY->SetMarkerStyle(7);
193     gptsXY->SetTitle(0);
194     gptsXY->Draw("AP");
195     gptsXY->GetXaxis()->SetTitle("Xglob");
196     gptsXY->GetYaxis()->SetTitle("Yglob");
197    }
198
199   TCanvas *cev1;
200   sprintf(ctit,"Event %d Zr",nevfordisp);
201   cev1=new TCanvas("cev1",ctit,600,600);
202   if(gptsRZ->GetN()>0){
203     gptsRZ->SetMarkerStyle(7);
204     gptsRZ->SetTitle(0);
205     gptsRZ->Draw("AP");
206     gptsRZ->GetXaxis()->SetTitle("Zglob");
207     gptsRZ->GetYaxis()->SetTitle("Radius");
208   }
209
210   return 0;
211 }
212
213