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 }