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 gAlice->GetRunLoader();
  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");
  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 gAlice->GetRunLoader();
 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");
 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 }