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