Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / ITS / ShowITSHitsRecPoints.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 ShowITSHitsRecPoints(Bool_t align=kFALSE,
20                            TString alignfile="ITSfullv11Misalignment.root")
21 {
22   ///////////////////////////////////////////////////////////////////////
23   // Macro to check clusters and hits in the 6 ITS layers              //
24   ///////////////////////////////////////////////////////////////////////
25
26   if (gClassTable->GetID("AliRun") < 0) {
27     gInterpreter->ExecuteMacro("loadlibs.C");
28   }
29   else { 
30     if(gAlice){
31       delete AliRunLoader::Instance();
32       delete gAlice;
33       gAlice=0;
34     }
35   }
36   // Set OCDB if needed
37   AliCDBManager* man = AliCDBManager::Instance();
38   if (!man->IsDefaultStorageSet()) {
39     printf("Setting a local default storage and run number 0\n");
40     man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
41     man->SetRun(0);
42   }
43   else {
44     printf("Using deafult storage \n");
45   }
46  
47   // retrives geometry 
48   if(!gGeoManager){
49     AliGeomManager::LoadGeometry("geometry.root");
50   }
51
52   if(align) {
53     TFile f(alignfile.Data());
54     TClonesArray* ar = (TClonesArray*)f.Get("ITSAlignObjs");
55     AliGeomManager::ApplyAlignObjsToGeom(*ar);
56     f.Close();
57   } else {
58     AliGeomManager::ApplyAlignObjsFromCDB("ITS");
59   }
60
61   AliRunLoader* rl = AliRunLoader::Open("galice.root");
62   if (rl == 0x0){
63     cerr<<"Can not open session RL=NULL"<< endl;
64     return -1;
65   }
66   Int_t retval = rl->LoadgAlice();
67   if (retval){
68     cerr<<"LoadgAlice returned error"<<endl;
69     return -1;
70   }
71   gAlice=rl->GetAliRun();
72
73   retval = rl->LoadHeader();
74   if (retval){
75     cerr<<"LoadHeader returned error"<<endl;
76     return -1;
77   }
78
79
80   AliITSLoader* ITSloader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
81   if(!ITSloader){
82     cerr<<"ITS loader not found"<<endl;
83     return -1;
84   }
85   ITSloader->LoadRecPoints("read");
86   ITSloader->LoadHits("read");
87
88
89   Float_t cluglo[3]={0.,0.,0.}; 
90   AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
91   ITS->SetTreeAddress();
92   AliITSgeom *geom = ITS->GetITSgeom();
93   AliITSDetTypeRec* detTypeRec = new AliITSDetTypeRec();
94   detTypeRec->SetITSgeom(ITSloader->GetITSgeom());
95   detTypeRec->SetDefaults();
96
97   Int_t modmin=geom->GetStartDet(0);
98   Int_t modmax=geom->GetLastDet(2);
99   Int_t totmod=modmax-modmin;
100   Float_t xlim[6]={4.5,7.5,16.,26.,40.,45.};
101   Float_t zlim[6]={15.,15.,22.,30.,45.,55.};
102
103   TH1F* hlayer=new TH1F("hlayer","",6,0.5,6.5);
104   TH1F** hmod=new TH1F*[6];
105   TH1F** hxl=new TH1F*[6];
106   TH1F** hzl=new TH1F*[6];
107   TH1F** hxg=new TH1F*[6];
108   TH1F** hyg=new TH1F*[6];
109   TH2F** hxyg=new TH2F*[6];
110   TH2F** hxygHits=new TH2F*[6];
111   TH1F** hzg=new TH1F*[6];
112   TH1F** hr=new TH1F*[6];
113   TH1F** hphi=new TH1F*[6];
114   TH1F** hq=new TH1F*[6];
115
116   Char_t name[10];
117   for(Int_t iLay=0;iLay<6;iLay++){
118     sprintf(name,"hmod%d",iLay+1);
119     hmod[iLay]=new TH1F(name,"",totmod,modmin-0.5,modmax+0.5);
120     hmod[iLay]->GetXaxis()->SetTitle("Module");
121     hmod[iLay]->GetXaxis()->CenterTitle();
122     sprintf(name,"hxloc%d",iLay+1);
123     hxl[iLay]=new TH1F(name,"",100,-4.,4.);
124     hxl[iLay]->GetXaxis()->SetTitle("Xloc (cm)");
125     hxl[iLay]->GetXaxis()->CenterTitle();
126     sprintf(name,"hzloc%d",iLay+1);
127     hzl[iLay]=new TH1F(name,"",100,-4.,4.);
128     hzl[iLay]->GetXaxis()->SetTitle("Zloc (cm)");
129     hzl[iLay]->GetXaxis()->CenterTitle();
130     sprintf(name,"hxgl%d",iLay+1);
131     hxg[iLay]=new TH1F(name,"",100,-xlim[iLay],xlim[iLay]);
132     hxg[iLay]->GetXaxis()->SetTitle("Xglob (cm)");
133     hxg[iLay]->GetXaxis()->CenterTitle();
134     sprintf(name,"hygl%d",iLay+1);
135     hyg[iLay]=new TH1F(name,"",100,-xlim[iLay],xlim[iLay]);
136     hyg[iLay]->GetXaxis()->SetTitle("Yglob (cm)");
137     hyg[iLay]->GetXaxis()->CenterTitle();
138     sprintf(name,"hxygl%d",iLay+1);
139     hxyg[iLay]=new TH2F(name,"",1000,-xlim[iLay],xlim[iLay],1000,-xlim[iLay],xlim[iLay]);
140     hxyg[iLay]->GetXaxis()->SetTitle("Xglob (cm)");
141     hxyg[iLay]->GetYaxis()->SetTitle("Yglob (cm)");
142     sprintf(name,"hxygHitsl%d",iLay+1);
143     hxygHits[iLay]=new TH2F(name,"",1000,-xlim[iLay],xlim[iLay],1000,-xlim[iLay],xlim[iLay]);
144     sprintf(name,"hzgl%d",iLay+1);
145     hzg[iLay]=new TH1F(name,"",100,-zlim[iLay],zlim[iLay]);
146     hzg[iLay]->GetXaxis()->SetTitle("Zglob (cm)");
147     hzg[iLay]->GetXaxis()->CenterTitle();
148     sprintf(name,"hr%d",iLay+1);
149     hr[iLay]=new TH1F(name,"",100,0.,50.);
150     hr[iLay]->GetXaxis()->SetTitle("r (cm)");
151     hr[iLay]->GetXaxis()->CenterTitle();
152     sprintf(name,"hphi%d",iLay+1);
153     hphi[iLay]=new TH1F(name,"",100,-TMath::Pi(),TMath::Pi());    
154     hphi[iLay]->GetXaxis()->SetTitle("#varphi (rad)");
155     hphi[iLay]->GetXaxis()->CenterTitle();
156     sprintf(name,"hq%d",iLay+1);
157     hq[iLay]=new TH1F(name,"",100,0.,300.);    
158     hq[iLay]->GetXaxis()->SetTitle("Charge (keV)");
159     hq[iLay]->GetXaxis()->CenterTitle();
160   }
161
162   TGraph **gpts=new TGraph*[3];
163   TGraph2D **gpts3d=new TGraph2D*[3];
164   TGraph **gRP=new TGraph*[6];
165   TGraph **gHend=new TGraph*[6];
166   TGraph **gHstart=new TGraph*[6];
167   for(Int_t i=0;i<3;i++){
168     gpts[i]=new TGraph(0);
169     gpts3d[i]=new TGraph2D(0);
170   }
171   for(Int_t i=0;i<6;i++){
172     gRP[i]=new TGraph(0);
173     gRP[i]->GetXaxis()->SetTitle("global x [cm]");
174     gRP[i]->GetYaxis()->SetTitle("global y [cm]");
175     gHend[i]=new TGraph(0);
176     gHstart[i]=new TGraph(0);
177   }
178   Int_t totev=rl->GetNumberOfEvents();
179   printf("Total Number of events = %d\n",totev);
180
181   Int_t iRP[6]={0},iH[6]={0};
182
183   for(Int_t iev=0;iev<totev;iev++){
184     rl->GetEvent(iev);
185     TTree *TR = ITSloader->TreeR();
186     TClonesArray *ITSrec  = detTypeRec->RecPoints();
187     TBranch *branch = 0;
188     if(TR && ITSrec){
189       branch = ITSloader->TreeR()->GetBranch("ITSRecPoints");
190       if(branch)branch->SetAddress(&ITSrec);
191     }
192     TTree *TH = ITSloader->TreeH();
193     TClonesArray *hits=new TClonesArray("AliITShit",10000);
194     TH->SetBranchAddress("ITS",&hits);
195
196     Int_t nparticles = rl->GetHeader()->GetNtrack();
197     cout<<"Event #"<<iev<<"   #Particles="<<nparticles<<endl;
198
199
200     Int_t ipt=0;
201     for(Int_t subd=0;subd<3;subd++){
202
203       Int_t first = geom->GetStartDet(subd);
204       Int_t last = geom->GetLastDet(subd);
205
206       for (Int_t mod=first; mod<=last; mod++){
207         detTypeRec->ResetRecPoints();
208         branch->GetEvent(mod);
209         Int_t nrecp = ITSrec->GetEntries();
210         if(nrecp>0){
211           for(Int_t irec=0;irec<nrecp;irec++) {
212             AliITSRecPoint *recp = (AliITSRecPoint*)ITSrec->At(irec);
213             Int_t lay=recp->GetLayer();
214             hlayer->Fill(lay);
215             recp->GetGlobalXYZ(cluglo);
216             Float_t rad=TMath::Sqrt(cluglo[0]*cluglo[0]+cluglo[1]*cluglo[1]); 
217             Float_t phi=TMath::ATan2(cluglo[1],cluglo[0]);
218             if(iev<3){
219               gpts[iev]->SetPoint(ipt,cluglo[0],cluglo[1]);
220               gpts3d[iev]->SetPoint(ipt,cluglo[0],cluglo[1],cluglo[2]);
221               ipt++;
222             }
223             gRP[lay]->SetPoint(iRP[lay],cluglo[0],cluglo[1]); iRP[lay]++; 
224             hmod[lay]->Fill(mod);
225             hzl[lay]->Fill(recp->GetDetLocalZ());
226             hxl[lay]->Fill(recp->GetDetLocalX());
227             hzg[lay]->Fill(cluglo[2]);
228             hyg[lay]->Fill(cluglo[1]);
229             hxyg[lay]->Fill(cluglo[0],cluglo[1]);
230             hxg[lay]->Fill(cluglo[0]);
231             hr[lay]->Fill(rad);
232             hphi[lay]->Fill(phi);
233             hq[lay]->Fill(recp->GetQ());
234           }
235         }
236       }
237     }
238
239
240     cout<<" Now read hits "<<endl;
241
242     Int_t nentrHits=(Int_t)TH->GetEntries();
243     for (Int_t i=0; i<nentrHits; i++) {
244       TH->GetEvent(i);
245       Int_t nhit=hits->GetEntriesFast();
246       for (Int_t ih=0; ih<nhit; ih++) {
247         AliITShit *h=(AliITShit*)hits->UncheckedAt(ih);
248         if(h->StatusExiting()) {
249           Double_t xl,yl,zl,tl,xl0,yl0,zl0,tl0;
250           h->GetPositionL(xl,yl,zl,tl);
251           h->GetPositionL0(xl0,yl0,zl0,tl0);
252           //if(TMath::Abs(yl-yl0)<0.0290) continue;
253           hxygHits[h->GetLayer()-1]->Fill(h->GetXG(),h->GetYG());
254           gHend[h->GetLayer()-1]->SetPoint(iH[h->GetLayer()-1],h->GetXG(),h->GetYG()); 
255           Double_t x0,y0,z0,t0;
256           h->GetPositionG0(x0,y0,z0,t0);
257           gHstart[h->GetLayer()-1]->SetPoint(iH[h->GetLayer()-1],x0,y0); 
258           iH[h->GetLayer()-1]++;
259           //printf("layer %d hit length xy %f hit length yl %f\n",h->GetLayer()-1,TMath::Sqrt((x0-h->GetXG())*(x0-h->GetXG())+(y0-h->GetYG())*(y0-h->GetYG())),TMath::Abs(yl-yl0));
260         } 
261       }
262     }
263   }
264
265
266   TCanvas **c=new TCanvas*[6];
267   Char_t ctit[30];
268   for(Int_t iLay=0;iLay<6;iLay++){
269     sprintf(name,"can%d",iLay+1);
270     sprintf(ctit,"Layer %d",iLay+1);
271     c[iLay]=new TCanvas(name,ctit,1200,900);
272     c[iLay]->Divide(3,3,0.001,0.001);
273     c[iLay]->cd(1);
274     hmod[iLay]->Draw();
275     c[iLay]->cd(2);
276     hxl[iLay]->Draw();
277     c[iLay]->cd(3);
278     hzl[iLay]->Draw();
279     c[iLay]->cd(4);
280     hxg[iLay]->Draw();
281     c[iLay]->cd(5);
282     hyg[iLay]->Draw();
283     c[iLay]->cd(6);
284     hzg[iLay]->Draw();
285     c[iLay]->cd(7);
286     hr[iLay]->Draw();
287     c[iLay]->cd(8);
288     hphi[iLay]->Draw();    
289     c[iLay]->cd(9);
290     hxyg[iLay]->Draw();    
291   }
292
293   TCanvas *cev0;
294   cev0=new TCanvas("cev0","Event 0",600,600);
295   gpts[0]->SetMarkerStyle(7);
296   gpts[0]->SetTitle(0);
297   gpts[0]->Draw("AP");
298
299   TCanvas *cev1;
300   cev1=new TCanvas("cev1","Event 1",600,600);
301   gpts[1]->SetMarkerStyle(7);
302   gpts[1]->SetTitle(0);
303   gpts[1]->Draw("AP");
304
305   TCanvas *cev2;
306   cev2=new TCanvas("cev2","Event 2",600,600);
307   gpts[2]->SetMarkerStyle(7);
308   gpts[2]->SetTitle(0);
309   gpts[2]->Draw("AP");
310
311   TCanvas *chr = new TCanvas("chr","chr");
312   chr->Divide(3,2);
313   for(Int_t i=0;i<6;i++) {
314     chr->cd(i+1);
315     gHend[i]->SetMarkerStyle(7);
316     gHend[i]->SetMarkerColor(1);
317     gHend[i]->Draw("A,P");
318     gHstart[i]->SetMarkerStyle(7);
319     gHstart[i]->SetMarkerColor(4);
320     gHstart[i]->Draw("P");
321     gRP[i]->SetMarkerStyle(7);
322     gRP[i]->SetMarkerColor(2);
323     gRP[i]->Draw("P");
324   }
325
326
327   return 0;
328 }
329
330
331 Int_t ShowITSHitsRecPointsNtuple(Bool_t align=kFALSE,
332                                  TString alignfile="ITSfullv11Misalignment.root")
333 {
334   ///////////////////////////////////////////////////////////////////////
335   // Macro to check clusters and hits in the 6 ITS layers              //
336   // Creates also ntuple                                               //
337   ///////////////////////////////////////////////////////////////////////
338
339   if (gClassTable->GetID("AliRun") < 0) {
340     gInterpreter->ExecuteMacro("loadlibs.C");
341   }
342   else { 
343     if(gAlice){
344       delete AliRunLoader::Instance();
345       delete gAlice;
346       gAlice=0;
347     }
348   }
349   // Set OCDB if needed
350   AliCDBManager* man = AliCDBManager::Instance();
351   if (!man->IsDefaultStorageSet()) {
352     printf("Setting a local default storage and run number 0\n");
353     man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
354     man->SetRun(0);
355   }
356   else {
357     printf("Using deafult storage \n");
358   }
359  
360   // retrives geometry 
361   if(!gGeoManager){
362     AliGeomManager::LoadGeometry("geometry.root");
363   }
364   if(align) {
365     TFile f(alignfile.Data());
366     TClonesArray* ar = (TClonesArray*)f.Get("ITSAlignObjs");
367     AliGeomManager::ApplyAlignObjsToGeom(*ar);
368     f.Close();
369   } else {
370     AliGeomManager::ApplyAlignObjsFromCDB("ITS");
371   }
372
373
374   AliRunLoader* rl = AliRunLoader::Open("galice.root");
375   if (rl == 0x0){
376     cerr<<"Can not open session RL=NULL"<< endl;
377     return -1;
378   }
379   Int_t retval = rl->LoadgAlice();
380   if (retval){
381     cerr<<"LoadgAlice returned error"<<endl;
382     return -1;
383   }
384   gAlice=rl->GetAliRun();
385
386   retval = rl->LoadHeader();
387   if (retval){
388     cerr<<"LoadHeader returned error"<<endl;
389     return -1;
390   }
391
392
393   AliITSLoader* ITSloader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
394   if(!ITSloader){
395     cerr<<"ITS loader not found"<<endl;
396     return -1;
397   }
398   ITSloader->LoadRecPoints("read");
399   ITSloader->LoadHits("read");
400
401
402   Float_t cluglo[3]={0.,0.,0.}; 
403   AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
404   ITS->SetTreeAddress();
405   AliITSgeom *geom = ITS->GetITSgeom();
406   AliITSDetTypeRec* detTypeRec = new AliITSDetTypeRec();
407   detTypeRec->SetITSgeom(ITSloader->GetITSgeom());
408   detTypeRec->SetDefaults();
409
410   Int_t modmin=geom->GetStartDet(0);
411   Int_t modmax=geom->GetLastDet(2);
412   Int_t totmod=modmax-modmin;
413   Float_t xlim[6]={4.5,7.5,16.,26.,40.,45.};
414   Float_t zlim[6]={15.,15.,22.,30.,45.,55.};
415
416   TH1F* hlayer=new TH1F("hlayer","",6,0.5,6.5);
417   TH1F** hmod=new TH1F*[6];
418   TH1F** hxl=new TH1F*[6];
419   TH1F** hzl=new TH1F*[6];
420   TH1F** hxg=new TH1F*[6];
421   TH1F** hyg=new TH1F*[6];
422   TH2F** hxyg=new TH2F*[6];
423   TH2F** hxygHits=new TH2F*[6];
424   TH1F** hzg=new TH1F*[6];
425   TH1F** hr=new TH1F*[6];
426   TH1F** hphi=new TH1F*[6];
427   TH1F** hq=new TH1F*[6];
428
429   Char_t name[10];
430   for(Int_t iLay=0;iLay<6;iLay++){
431     sprintf(name,"hmod%d",iLay+1);
432     hmod[iLay]=new TH1F(name,"",totmod,modmin-0.5,modmax+0.5);
433     hmod[iLay]->GetXaxis()->SetTitle("Module");
434     hmod[iLay]->GetXaxis()->CenterTitle();
435     sprintf(name,"hxloc%d",iLay+1);
436     hxl[iLay]=new TH1F(name,"",100,-4.,4.);
437     hxl[iLay]->GetXaxis()->SetTitle("Xloc (cm)");
438     hxl[iLay]->GetXaxis()->CenterTitle();
439     sprintf(name,"hzloc%d",iLay+1);
440     hzl[iLay]=new TH1F(name,"",100,-4.,4.);
441     hzl[iLay]->GetXaxis()->SetTitle("Zloc (cm)");
442     hzl[iLay]->GetXaxis()->CenterTitle();
443     sprintf(name,"hxgl%d",iLay+1);
444     hxg[iLay]=new TH1F(name,"",100,-xlim[iLay],xlim[iLay]);
445     hxg[iLay]->GetXaxis()->SetTitle("Xglob (cm)");
446     hxg[iLay]->GetXaxis()->CenterTitle();
447     sprintf(name,"hygl%d",iLay+1);
448     hyg[iLay]=new TH1F(name,"",100,-xlim[iLay],xlim[iLay]);
449     hyg[iLay]->GetXaxis()->SetTitle("Yglob (cm)");
450     hyg[iLay]->GetXaxis()->CenterTitle();
451     sprintf(name,"hxygl%d",iLay+1);
452     hxyg[iLay]=new TH2F(name,"",1000,-xlim[iLay],xlim[iLay],1000,-xlim[iLay],xlim[iLay]);
453     hxyg[iLay]->GetXaxis()->SetTitle("Xglob (cm)");
454     hxyg[iLay]->GetYaxis()->SetTitle("Yglob (cm)");
455     sprintf(name,"hxygHitsl%d",iLay+1);
456     hxygHits[iLay]=new TH2F(name,"",1000,-xlim[iLay],xlim[iLay],1000,-xlim[iLay],xlim[iLay]);
457     sprintf(name,"hzgl%d",iLay+1);
458     hzg[iLay]=new TH1F(name,"",100,-zlim[iLay],zlim[iLay]);
459     hzg[iLay]->GetXaxis()->SetTitle("Zglob (cm)");
460     hzg[iLay]->GetXaxis()->CenterTitle();
461     sprintf(name,"hr%d",iLay+1);
462     hr[iLay]=new TH1F(name,"",100,0.,50.);
463     hr[iLay]->GetXaxis()->SetTitle("r (cm)");
464     hr[iLay]->GetXaxis()->CenterTitle();
465     sprintf(name,"hphi%d",iLay+1);
466     hphi[iLay]=new TH1F(name,"",100,-TMath::Pi(),TMath::Pi());    
467     hphi[iLay]->GetXaxis()->SetTitle("#varphi (rad)");
468     hphi[iLay]->GetXaxis()->CenterTitle();
469     sprintf(name,"hq%d",iLay+1);
470     hq[iLay]=new TH1F(name,"",100,0.,300.);    
471     hq[iLay]->GetXaxis()->SetTitle("Charge (keV)");
472     hq[iLay]->GetXaxis()->CenterTitle();
473   }
474
475   TGraph **gpts=new TGraph*[3];
476   TGraph2D **gpts3d=new TGraph2D*[3];
477   TGraph **gRP=new TGraph*[6];
478   TGraph **gHend=new TGraph*[6];
479   TGraph **gHstart=new TGraph*[6];
480   for(Int_t i=0;i<3;i++){
481     gpts[i]=new TGraph(0);
482     gpts3d[i]=new TGraph2D(0);
483   }
484   for(Int_t i=0;i<6;i++){
485     gRP[i]=new TGraph(0);
486     gRP[i]->GetXaxis()->SetTitle("global x [cm]");
487     gRP[i]->GetYaxis()->SetTitle("global y [cm]");
488     gHend[i]=new TGraph(0);
489     gHstart[i]=new TGraph(0);
490   }
491   Int_t totev=rl->GetNumberOfEvents();
492   printf("Total Number of events = %d\n",totev);
493
494   Int_t iRP[6]={0},iH[6]={0};
495
496   TNtuple *nt = new TNtuple("nt","ntuple","lay:mod:deltaxl:deltax:deltay:deltaz:phi:z:xl:hlength");
497
498   for(Int_t iev=0;iev<totev;iev++){
499     rl->GetEvent(iev);
500     TTree *TR = ITSloader->TreeR();
501     TClonesArray *ITSrec  = detTypeRec->RecPoints();
502     TBranch *branch = 0;
503     if(TR && ITSrec){
504       branch = ITSloader->TreeR()->GetBranch("ITSRecPoints");
505       if(branch)branch->SetAddress(&ITSrec);
506     }
507     TTree *TH = ITSloader->TreeH();
508     TClonesArray *hits=new TClonesArray("AliITShit",10000);
509     TH->SetBranchAddress("ITS",&hits);
510
511     Int_t nparticles = rl->GetHeader()->GetNtrack();
512     cout<<"Event #"<<iev<<"   #Particles="<<nparticles<<endl;
513
514
515     Int_t ipt=0;
516     for(Int_t subd=0;subd<3;subd++){
517
518       Int_t first = geom->GetStartDet(subd);
519       Int_t last = geom->GetLastDet(subd);
520
521       for (Int_t mod=first; mod<=last; mod++){
522         detTypeRec->ResetRecPoints();
523         branch->GetEvent(mod);
524         Int_t nrecp = ITSrec->GetEntries();
525         if(nrecp>0){
526           for(Int_t irec=0;irec<nrecp;irec++) {
527             AliITSRecPoint *recp = (AliITSRecPoint*)ITSrec->At(irec);
528             Int_t lay=recp->GetLayer();
529             hlayer->Fill(lay);
530             recp->GetGlobalXYZ(cluglo);
531             Float_t rad=TMath::Sqrt(cluglo[0]*cluglo[0]+cluglo[1]*cluglo[1]); 
532             Float_t phi=TMath::ATan2(cluglo[1],cluglo[0]);
533             if(cluglo[2]>0) {gRP[lay]->SetPoint(iRP[lay],cluglo[0],cluglo[1]); iRP[lay]++;} 
534             hmod[lay]->Fill(mod);
535             hzl[lay]->Fill(recp->GetDetLocalZ());
536             hxl[lay]->Fill(recp->GetDetLocalX());
537             hzg[lay]->Fill(cluglo[2]);
538             hyg[lay]->Fill(cluglo[1]);
539             hxyg[lay]->Fill(cluglo[0],cluglo[1]);
540             hxg[lay]->Fill(cluglo[0]);
541             hr[lay]->Fill(rad);
542             hphi[lay]->Fill(phi);
543             hq[lay]->Fill(recp->GetQ());
544
545             Double_t hlength=0.;
546             for(Int_t jhits=0;jhits<TH->GetEntries();jhits++) {
547               TH->GetEvent(jhits);
548               Int_t nhit=hits->GetEntriesFast();
549               for (Int_t ih=0; ih<nhit; ih++) {
550                 AliITShit *h=(AliITShit*)hits->UncheckedAt(ih);
551                 if(h->GetTrack()!=recp->GetLabel(0)) continue;
552                 if(h->GetLayer()-1!=lay) continue;
553                 if(h->GetModule()!=mod) continue;
554                 //hxygHits[h->GetLayer()-1]->Fill(h->GetXG(),h->GetYG());
555                 //gHend[h->GetLayer()-1]->SetPoint(iH[h->GetLayer()-1],h->GetXG(),h->GetYG()); 
556                 Double_t x0,y0,z0,t0;
557                 Double_t x,y,z,t;
558                 h->GetPositionG0(x0,y0,z0,t0);
559                 h->GetPositionG(x,y,z,t);
560                 //gHstart[h->GetLayer()-1]->SetPoint(iH[h->GetLayer()-1],x0,y0); 
561                 //iH[h->GetLayer()-1]++;
562                 Double_t xl,yl,zl,tl,xl0,yl0,zl0,tl0;
563                 h->GetPositionL(xl,yl,zl,tl);
564                 h->GetPositionL0(xl0,yl0,zl0,tl0);
565                 hlength += TMath::Abs(yl-yl0);
566                 if(!h->StatusExiting()) continue;
567                 nt->Fill(lay,
568                          mod,
569                          1.e4*(0.5*(xl+xl0)-recp->GetDetLocalX()),
570                          1.e4*(0.5*(x+x0)-cluglo[0]),1.e4*(0.5*(y+y0)-cluglo[1]),
571                          1.e4*(0.5*(z+z0)-cluglo[2]),
572                          phi,
573                          cluglo[2],
574                          recp->GetDetLocalX(),
575                          hlength);
576                 hlength = 0.;
577               }
578             }
579
580
581           }
582         }
583       }
584     }
585
586
587   }
588
589
590   TCanvas **c=new TCanvas*[6];
591   Char_t ctit[30];
592   for(Int_t iLay=0;iLay<6;iLay++){
593     sprintf(name,"can%d",iLay+1);
594     sprintf(ctit,"Layer %d",iLay+1);
595     c[iLay]=new TCanvas(name,ctit,1200,900);
596     c[iLay]->Divide(3,3,0.001,0.001);
597     c[iLay]->cd(1);
598     hmod[iLay]->Draw();
599     c[iLay]->cd(2);
600     hxl[iLay]->Draw();
601     c[iLay]->cd(3);
602     hzl[iLay]->Draw();
603     c[iLay]->cd(4);
604     hxg[iLay]->Draw();
605     c[iLay]->cd(5);
606     hyg[iLay]->Draw();
607     c[iLay]->cd(6);
608     hzg[iLay]->Draw();
609     c[iLay]->cd(7);
610     hr[iLay]->Draw();
611     c[iLay]->cd(8);
612     hphi[iLay]->Draw();    
613     c[iLay]->cd(9);
614     hxyg[iLay]->Draw();    
615   }
616
617   TCanvas *cev0;
618   cev0=new TCanvas("cev0","Event 0",600,600);
619   gpts[0]->SetMarkerStyle(7);
620   gpts[0]->SetTitle(0);
621   gpts[0]->Draw("AP");
622
623   TCanvas *cev1;
624   cev1=new TCanvas("cev1","Event 1",600,600);
625   gpts[1]->SetMarkerStyle(7);
626   gpts[1]->SetTitle(0);
627   gpts[1]->Draw("AP");
628
629   TCanvas *cev2;
630   cev2=new TCanvas("cev2","Event 2",600,600);
631   gpts[2]->SetMarkerStyle(7);
632   gpts[2]->SetTitle(0);
633   gpts[2]->Draw("AP");
634
635   TFile *out= new TFile("ntHitsRecPoints.root","recreate");
636   nt->Write();
637   out->Close();
638
639   return 0;
640 }