ITS/AliITSGeoPlot.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include<Riostream.h>
   3 #include<TROOT.h>
   4 #include<TArrayI.h>
   5 #include<TBranch.h>
   6 #include<TCanvas.h>
   7 #include<TClassTable.h>
   8 #include<TClonesArray.h>
   9 #include<TFile.h>
  10 #include<TH1.h>
  11 #include<TH2.h>
  12 #include<TObject.h>
  13 #include<TObjArray.h>
  14 #include<TTree.h>
  15 #include "AliGenEventHeader.h"
  16 #include <AliRun.h>
  17 #include <AliITS.h>
  18 #include <AliITSgeom.h>
  19 #include <AliITSDetType.h>
  20 #include <AliITSRecPoint.h>
  21 #include <AliITSdigit.h>
  22 #include <AliITShit.h>
  23 #include <AliITSmodule.h>
  24 #include <AliITSsegmentation.h>
  25 #include <AliITSsegmentationSPD.h>
  26 #include <AliITSsegmentationSDD.h>
  27 #include <AliITSsegmentationSSD.h>
  28 #include <AliRunLoader.h>
  29 #include <AliITSLoader.h>
  30 #include <AliHeader.h>
  31 #endif
  32 void GetHitsCoor(TObject *its, Int_t mod, TObjArray & histos, Int_t subd,Bool_t verb);
  33 Int_t GetRecCoor(TObject *ge, TClonesArray *ITSrec, Int_t mod, TH2F *h2, TH1F *h1, Bool_t verb);
  34 void GetDigits(TObject *tmps,TObject *ge,TClonesArray *ITSdigits, Int_t subd, Int_t mod, Bool_t verbose, TObjArray & histos);
  35
  36 Int_t AliITSGeoPlot (Int_t evesel=0, char *opt="All+Rec", char *filename="galice.root", Int_t isfastpoints = 0) {
  37   /*******************************************************************
  38    *  This macro displays geometrical information related to the
  39    *  hits, digits and rec points in ITS.
  40    *  There are histograms that are not displayed (i.e. energy
  41    *  deposition) but that are saved onto a file (see below)
  42    *
  43    *  INPUT arguments:
  44    *
  45    *  Int_t evesel:  the selected event number. It's =0 by default
  46    *
  47    *  Options: Any combination of:
  48    *    1) subdetector name:  "SPD", "SDD", "SSD", "All" (default)
  49    *    2) Printouts:        "Verbose" Almost everything is printed out:
  50    *                          it is wise to redirect the output onto a file
  51    *                    e.g.: .x AliITSGeoPlot.C("All+Verbose") > out.log
  52    *    3) Rec Points option: "Rec"   ---> Uses Rec. Points (default)
  53    *                           otherwise ---> uses hits and digits only
  54    *    Examples:
  55    *       .x AliITSGeoPlot();  (All subdetectors; no-verbose; no-recpoints)
  56    *       .x AliITSGeoPlot("SPD+SSD+Verbose+Rec");
  57    *
  58    *    filename:   It's "galice.root" by default.
  59    *    isfastpoints: integer. It is set to 0 by defaults. This means that
  60    *                slow reconstruction is assumed. If fast recpoint are in
  61    *                in use, isfastpoints must be set =1.
  62    *    BEWARE:     to analyze digits/recpoints  produced with aliroot
  63    *                versions older than v3-07-Release, the input parameters
  64    *                filename, FileDigits, FileRec and isfastpoints must be left
  65    *                to their default value
  66    *
  67    *  OUTPUT: It produces a root file with a list of histograms
  68    *
  69    *  WARNING: spatial information for SSD/DIGITS is obtained by pairing
  70    *          digits on p and n side originating from the same track, when
  71    *          possible. This (mis)use of DIGITS is tolerated for debugging
  72    *          purposes only !!!!  The pairing in real life should be done
  73    *          starting from info really available...
  74    *
  75    *  COMPILATION: this macro can be compiled.
  76    *      1)       You need to set your include path with
  77    * gSystem->SetIncludePath("-I- -I$ALICE_ROOT/include -I$ALICE_ROOT/ITS -g");
  78    *      3)       If you are using root instead of aliroot you need to
  79    *               execute the macro loadlibs.C in advance
  80    *      4)       To compile this macro from root (or aliroot):
  81    *                 ---  .L AliITSGeoPlot.C++
  82    *                 ---  AliITSGeoPlot();
  83    *
  84    *  M.Masera  14/05/2001 18:30
  85    *  Last rev. 09/06/2003 17:00 (Adapted to NewIO)  m.m.
  86    ********************************************************************/
  87
  88   //Options
  89   TString choice(opt);
  90   Bool_t All = choice.Contains("All");
  91   Bool_t verbose=choice.Contains("Verbose");
  92   Bool_t userec=choice.Contains("Rec");
  93   Int_t retcode=1; //return code
  94 #if !(!defined(__CINT__) || defined(__MAKECINT__))
  95   if (gClassTable->GetID("AliRun") < 0) {
  96     gROOT->LoadMacro("loadlibs.C");
  97     loadlibs();
  98   }
  99   else {
 100 #endif
 101     if(gAlice){
 102       delete gAlice->GetRunLoader();
 103       delete gAlice;
 104       gAlice=0;
 105     }
 106 #if !(!defined(__CINT__) || defined(__MAKECINT__))
 107   }
 108 #endif
 109
 110 AliRunLoader* rl = AliRunLoader::Open(filename);
 111  if (rl == 0x0){
 112    cerr<<"AliITSGeoPlot.C : Can not open session RL=NULL"<< endl;
 113    return -1;
 114  }
 115  Int_t retval = rl->LoadgAlice();
 116  if (retval){
 117    cerr<<"AliITSGeoPlot.C : LoadgAlice returned error"<<endl;
 118    return -1;
 119  }
 120  gAlice=rl->GetAliRun();
 121
 122  retval = rl->LoadHeader();
 123  if (retval){
 124    cerr<<"AliITSGeoPlot.C : LoadHeader returned error"<<endl;
 125    return -1;
 126  }
 127
 128  AliITSLoader* ITSloader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
 129
 130  if(!ITSloader){
 131    cerr<<"AliITSGeoPlot.C :  ITS loader not found"<<endl;
 132    return -1;
 133  }
 134
 135  ITSloader->LoadHits("read");
 136  ITSloader->LoadDigits("read");
 137  if(isfastpoints==1)ITSloader->SetRecPointsFileName("ITS.FastRecPoints.root");
 138  ITSloader->LoadRecPoints("read");
 139  rl->GetEvent(evesel);
 140  Int_t nparticles = rl->GetHeader()->GetNtrack();
 141  AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
 142  ITS->SetTreeAddress();
 143   if(verbose) {
 144     cout<<" "<<endl<<" "<<endl;
 145     cout<<"******* Event processing started   *******"<<endl;
 146     cout<<"In the following, hits with 'StatusEntering' flag active"<<endl;
 147     cout<<"will not be processed"<<endl;
 148     cout << "Number of particles=  " << nparticles <<endl;
 149   }
 150
 151   // HITS
 152   TTree *TH = ITSloader->TreeH();
 153   Stat_t ntracks = TH->GetEntries();
 154   if(verbose)cout<<"Number of primary tracks= "<<ntracks<<endl;
 155
 156   // ITS
 157   Int_t nmodules;
 158   ITS->InitModules(-1,nmodules);
 159   cout<<"Number of ITS modules= "<<nmodules<<endl;
 160   cout<<"Filling modules... It takes a while, now. Please be patient"<<endl;
 161   ITS->FillModules(0,0,nmodules," "," ");
 162   cout<<"ITS modules .... DONE!"<<endl;
 163
 164   // DIGITS
 165   TTree *TD = ITSloader->TreeD();
 166
 167   //RECPOINTS
 168   TTree *TR = ITSloader->TreeR();
 169   TClonesArray *ITSrec  = ITS->RecPoints();
 170   TBranch *branch = 0;
 171   if(userec && TR && ITSrec){
 172     if(isfastpoints==1){
 173       branch = ITSloader->TreeR()->GetBranch("ITSRecPointsF");
 174       cout<<"using fast points\n";
 175     }
 176     else {
 177       branch = ITSloader->TreeR()->GetBranch("ITSRecPoints");
 178     }
 179     if(branch)branch->SetAddress(&ITSrec);
 180   }
 181
 182   if(userec && (!TR || !ITSrec || !branch)){
 183     userec = kFALSE;
 184     cout<<"\n ======================================================= \n";
 185     cout<<"WARNING: there are no RECPOINTS on this file ! \n";
 186     cout<<"======================================================= \n \n";
 187   }
 188
 189   //local variables
 190   Int_t mod;   //module number
 191   Int_t nbytes = 0;
 192   Double_t pos[3];  // Global position of the current module
 193   Float_t ragdet; // Radius of detector (x y plane)
 194   Int_t first; // first module
 195   Int_t last; // last module
 196   Int_t nrecp; //number of RecPoints for a given module
 197
 198   //List of histograms
 199   TObjArray histos(26,0);  // contains the pointers to the histograms
 200   // Book histograms SPD
 201   TH2F *rspd = new TH2F("rspd","Radii of digits - SPD",50,-10.,10.,50,-10.,10.);
 202   TH2F *rhspd = new TH2F("rhspd","Radii of hits - SPD",50,-10.,10.,50,-10.,10.);
 203   TH2F *rmspd = new TH2F("rmspd","Radii of SPD modules",50,-10.,10.,50,-10.,10.);
 204   TH1F *zspd = new TH1F("zspd","Z of digits - SPD",100,-30.,30.);
 205   TH1F *zhspd = new TH1F("zhspd","Z of hits - SPD",100,-30.,30.);
 206   TH1F *zmspd = new TH1F("zmspd","Z of SPD modules",100,-30,30.);
 207   TH2F *rrspd = new TH2F("rrspd","Radii of recpoints - SPD",50,-10.,10.,50,-10.,10.);
 208   TH1F *zrspd = new TH1F("zrspd","Z of recpoints - SPD",100,-30.,30.);
 209   TH1F *enespd = new TH1F("enespd","Energy deposition SPD (KeV)",100,0.,1000.);
 210   histos.AddLast(rspd);  // 0
 211   histos.AddLast(rhspd); // 1
 212   histos.AddLast(rmspd); // 2
 213   histos.AddLast(zspd);  // 3
 214   histos.AddLast(zhspd); // 4
 215   histos.AddLast(zmspd); // 5
 216   histos.AddLast(rrspd); // 6
 217   histos.AddLast(zrspd); // 7
 218   histos.AddLast(enespd); // 8
 219   // Book histograms SDD
 220   TH2F *rsdd = new TH2F("rsdd","Radii of digits - SDD",50,-40.,40.,50,-40.,40.);
 221   TH2F *rhsdd = new TH2F("rhsdd","Radii of hits - SDD",50,-40.,40.,50,-40.,40.);
 222   TH2F *rmsdd = new TH2F("rmsdd","Radii of SDD modules",50,-40.,40.,50,-40.,40.);
 223   TH1F *zsdd = new TH1F("zsdd","Z of digits - SDD",100,-40.,40.);
 224   TH1F *zhsdd = new TH1F("zhsdd","Z of hits - SDD",100,-40.,40.);
 225   TH1F *zmsdd = new TH1F("zmsdd","Z of SDD modules",100,-40,40.);
 226   TH2F *rrsdd = new TH2F("rrsdd","Radii of recpoints - SDD",50,-40.,40.,50,-40.,40.);
 227   TH1F *zrsdd = new TH1F("zrsdd","Z of recpoints - SDD",100,-40.,40.);
 228   TH1F *enesdd = new TH1F("enesdd","Energy deposition SDD (KeV)",100,0.,1000.);
 229   histos.AddLast(rsdd);  // 9
 230   histos.AddLast(rhsdd); // 10
 231   histos.AddLast(rmsdd); // 11
 232   histos.AddLast(zsdd);  // 12
 233   histos.AddLast(zhsdd); // 13
 234   histos.AddLast(zmsdd); // 14
 235   histos.AddLast(rrsdd); // 15
 236   histos.AddLast(zrsdd); // 16
 237   histos.AddLast(enesdd); // 17
 238   // Book histogram SSD
 239   TH2F *rssd = new TH2F("rssd","Radii of digits - SSD",50,-50.,50.,50,-50.,50.);
 240   TH2F *rhssd = new TH2F("rhssd","Radii of hits - SSD",50,-50.,50.,50,-50.,50.);
 241   TH2F *rmssd = new TH2F("rmssd","Radii of SSD modules",50,-50.,50.,50,-50.,50.);
 242   TH1F *zssd = new TH1F("zssd","Z of digits - SSD",100,-70.,70.);
 243   TH1F *zhssd = new TH1F("zhssd","Z of hits - SSD",100,-70.,70.);
 244   TH1F *zmssd = new TH1F("zmssd","Z of SSD modules",100,-70,70.);
 245   TH2F *rrssd = new TH2F("rrssd","Radii of recpoints - SSD",50,-50.,50.,50,-50.,50.);
 246   TH1F *zrssd = new TH1F("zrssd","Z of recpoints - SSD",100,-70.,70.);
 247   TH1F *enessd = new TH1F("enessd","Energy deposition SSD (KeV)",100,0.,1000.);
 248   histos.AddLast(rssd);  // 18
 249   histos.AddLast(rhssd); // 19
 250   histos.AddLast(rmssd); // 20
 251   histos.AddLast(zssd);  // 21
 252   histos.AddLast(zhssd); // 22
 253   histos.AddLast(zmssd); // 23
 254   histos.AddLast(rrssd); // 24
 255   histos.AddLast(zrssd); // 25
 256   histos.AddLast(enessd); // 26
 257   //
 258   // Loop on subdetectors
 259   //
 260   AliITSgeom *geom = ITS->GetITSgeom();
 261   TString detna; // subdetector name
 262   for(Int_t subd=0;subd<3;subd++){
 263     if(All || (choice.Contains("SPD") && subd==0) || (choice.Contains("SDD") && subd==1) || (choice.Contains("SSD") && subd==2)){
 264       // Prepare array for the digits
 265       TClonesArray *ITSdigits  = ITS->DigitsAddress(subd);
 266       Bool_t usedigits = kTRUE;
 267       if(!ITSdigits){
 268         usedigits = kFALSE;
 269         cout<<"\n ======================================================= \n";
 270         cout<<"WARNING: there are no DIGITS on this file ! \n";
 271         cout<<"======================================================= \n \n";
 272       }
 273       // Get segmentation model
 274       AliITSDetType *iDetType=ITS->DetType(subd);
 275       if(subd==0)detna="SPD";
 276       if(subd==1)detna="SDD";
 277       if(subd==2)detna="SSD";
 278       AliITSsegmentation *seg=(AliITSsegmentation*)iDetType->GetSegmentationModel();
 279       // Loop on modules
 280       first = geom->GetStartDet(subd);
 281       last = geom->GetLastDet(subd);
 282       if(verbose){
 283         cout<<"     "<<endl<<"-------------------------------------"<<endl;
 284         cout<<"Start processing subdetector "<<detna<<endl;
 285         cout<<detna<<" first module "<<first<<endl;
 286         cout<<detna<<" last module "<<last<<endl;
 287         cout<<" "<<endl<<" "<<endl;
 288       }
 289       for (mod=first; mod<=last; mod++){
 290         geom->GetTrans(mod,pos);  // position of the module in the MRS
 291         ragdet=sqrt(pos[0]*pos[0]+pos[1]*pos[1]);
 292         // The following 2 histos are a check of the geometry
 293         TH2F *bidi = (TH2F*)histos.At(2+subd*9);
 294         TH1F *uni = (TH1F*)histos.At(5+subd*9);
 295         bidi->Fill(pos[0],pos[1]);
 296         uni->Fill(pos[2]);
 297         if(verbose){
 298           cout<<"=========================================================\n";
 299           cout<<detna<<" module="<<mod<<endl;
 300           cout<<"Mod. coordinates: "<<pos[0]<<", "<<pos[1]<<", ";
 301           cout<<pos[2]<<" Radius "<<ragdet<<endl;
 302         }
 303
 304         // Hits
 305         GetHitsCoor(ITS,mod,histos,subd,verbose);
 306
 307         //RecPoints
 308         if(userec){
 309           ITS->ResetRecPoints();
 310           //          TR->GetEvent(mod);
 311           branch->GetEvent(mod);
 312           TH2F *bidi=(TH2F*)histos.At(6+subd*9);
 313           TH1F *uni=(TH1F*)histos.At(7+subd*9);
 314           nrecp=GetRecCoor(geom,ITSrec,mod,bidi,uni,verbose);
 315         }
 316
 317         // Digits
 318         if(usedigits){
 319           ITS->ResetDigits();
 320           nbytes += TD->GetEvent(mod);
 321           GetDigits(seg,geom,ITSdigits,subd,mod,verbose,histos);
 322         }
 323
 324       } // End of loop on the modules
 325       TH1F *h1tmp;
 326       TH2F *h2tmp;
 327       //  Plot the histograms
 328       TCanvas *current=0; // current Canvas (1--> SPD, 2---> SDD, 3---> SSD)
 329       if(subd==0){
 330         // Prepare canvas 1
 331         TCanvas *c1 = new TCanvas("c1","SPD",10,10,600,900);
 332         c1->Divide(2,3);
 333         current=c1;
 334       }
 335       if(subd==1){
 336         // Prepare canvas 2
 337         TCanvas *c2 = new TCanvas("c2","SDD",40,40,600,900);
 338         c2->Divide(2,3);
 339         current=c2;
 340       }
 341       if(subd==2){
 342         // Prepare canvas 3
 343         TCanvas *c3 = new TCanvas("c3","SSD",70,70,600,900);
 344         c3->Divide(2,3);
 345         current=c3;
 346       }
 347       current->cd(1);
 348       h2tmp = (TH2F*)histos.At(9*subd);
 349       h2tmp->Draw();
 350       current->cd(2);
 351       h1tmp=(TH1F*)histos.At(3+subd*9);
 352       h1tmp->Draw();
 353       current->cd(3);
 354       h2tmp=(TH2F*)histos.At(1+9*subd);
 355       h2tmp->Draw();
 356       current->cd(4);
 357       h1tmp=(TH1F*)histos.At(4+subd*9);
 358       h1tmp->Draw();
 359
 360       if(userec){
 361         current->cd(5);
 362         h2tmp=(TH2F*)histos.At(6+9*subd);
 363         h2tmp->Draw();
 364         current->cd(6);
 365         h1tmp=(TH1F*)histos.At(7+subd*9);
 366         h1tmp->Draw();
 367       }
 368
 369       else {
 370         current->cd(5);
 371         h2tmp=(TH2F*)histos.At(2+9*subd);
 372         h2tmp->Draw();
 373         current->cd(6);
 374         h2tmp=(TH2F*)histos.At(5+9*subd);
 375         h2tmp->Draw();
 376       }
 377     } // if(All.....
 378   } // end of loop on subdetectors
 379   // Save the histograms
 380   TFile *fh = new TFile("AliITSGeoPlot.root","recreate");
 381   // The list is written to file as a single entry
 382   TList *lihi = new TList();
 383   // copy the pointers to the histograms to a TList object.
 384   // The histograms concerning recpoints are not copied if
 385   // 'userec' is false.
 386   for(Int_t i=0;i<histos.GetEntriesFast();i++){
 387     if(choice.Contains("All") || (choice.Contains("SPD") && i<8) || (choice.Contains("SDD") && i>7 && i<16) || (choice.Contains("SSD") && i>15)){
 388       if(!(!userec && ((i+2)%9==0 || (i+1)%9==0)))lihi->Add(histos.At(i));
 389     }
 390   }
 391   lihi->Write("Histograms ITS hits+digits+recpoints",TObject::kSingleKey);
 392   fh->Close();
 393
 394   return retcode;
 395 }
 396
 397
 398 void GetHitsCoor(TObject *its, Int_t mod, TObjArray & histos, Int_t subd,Bool_t verb){
 399   TH2F *h2=(TH2F*)histos.At(1+subd*9);
 400   TH1F *h1=(TH1F*)histos.At(4+subd*9);
 401   TH1F *ener=(TH1F*)histos.At(8+subd*9);
 402   AliITS *ITS= (AliITS*)its;
 403   AliITSmodule *modu = ITS->GetModule(mod);
 404   TObjArray *fHits = modu->GetHits();
 405   Int_t nhits = fHits->GetEntriesFast();
 406   if(nhits>0){
 407     if(verb){
 408       cout<<"-------------------------------------------------------"<<endl;
 409       cout<<"Number of HITS for module "<<mod<<": "<<nhits<<endl;
 410     }
 411     for (Int_t hit=0;hit<nhits;hit++) {
 412       AliITShit *iHit = (AliITShit*) fHits->At(hit);
 413       if(!iHit->StatusEntering()){
 414         Float_t x=iHit->GetXG();
 415         Float_t y=iHit->GetYG();
 416         Float_t z=iHit->GetZG();
 417         Float_t edep=iHit->GetIonization()*1000000;
 418         h2->Fill(x,y);
 419         h1->Fill(z);
 420         ener->Fill(edep);
 421         if(verb){
 422           cout<<"hit # "<<hit<<" Global coordinates "<<x<<" "<<y<<" "<<z<<endl;
 423           cout<<"track # "<<iHit->GetTrack()<<" energy deposition (KeV)= ";
 424           cout<<edep<<endl;
 425         }
 426       }
 427     }
 428   }
 429 }
 430
 431
 432 Int_t GetRecCoor(TObject *ge, TClonesArray *ITSrec, Int_t mod, TH2F *h2, TH1F *h1, Bool_t verb){
 433   AliITSgeom *geom = (AliITSgeom*)ge;
 434   Int_t nrecp = ITSrec->GetEntries();
 435   if(nrecp>0){
 436     Float_t lc[3]; for(Int_t i=0; i<3; i++) lc[i]=0.;
 437     Float_t gc[3]; for(Int_t i=0; i<3; i++) gc[i]=0.;
 438     if(verb){
 439       cout<<"-------------------------------------------------------"<<endl;
 440       cout<<"Number of REC POINTS for module "<<mod<<": "<<nrecp<<endl;
 441     }
 442     for(Int_t irec=0;irec<nrecp;irec++) {
 443       AliITSRecPoint *recp = (AliITSRecPoint*)ITSrec->At(irec);
 444       lc[0]=recp->GetX();
 445       lc[2]=recp->GetZ();
 446       geom->LtoG(mod,lc,gc);
 447       if(verb){
 448         cout<<"recp # "<<irec<<" local coordinates. lx= "<<lc[0]<<" lz= ";
 449         cout<<lc[2]<<endl;
 450         Float_t r=sqrt(gc[0]*gc[0]+gc[1]*gc[1]);
 451         cout<<"Global coor. + radius "<<gc[0]<<" "<<gc[1]<<" "<<gc[2]<<" ";
 452         cout<<r<<endl;
 453         cout<<"Associated track "<<recp->GetLabel(0)<<endl;
 454       }
 455       h2->Fill(gc[0],gc[1]);
 456       h1->Fill(gc[2]);
 457     }
 458   }
 459   return nrecp;
 460 }
 461
 462 void GetDigits(TObject *tmps,TObject *ge,TClonesArray *ITSdigits, Int_t subd, Int_t mod, Bool_t verbose, TObjArray & histos){
 463   AliITSsegmentation *seg = (AliITSsegmentation*)tmps;
 464   AliITSgeom *geom = (AliITSgeom*)ge;
 465   Int_t layer;
 466   Int_t ladder;
 467   Int_t detec;
 468   if(subd==2){
 469     geom->GetModuleId(mod,layer,ladder,detec);
 470     seg->SetLayer(layer);
 471   }
 472   Float_t lcoor[3]; for(Int_t j=0; j<3; j++) lcoor[j]=0.;  //local coord dig.
 473   Float_t gcoor[3]; for(Int_t j=0; j<3; j++) gcoor[j]=0.; // global coo. dig.
 474   Float_t ragdig; // Radius digit
 475   TArrayI ssdone(5000);  // used to match p and n side digits of strips
 476   TArrayI pair(5000);    // as above
 477   Int_t ndigits = ITSdigits->GetEntries();
 478   AliITSdigit *digs;
 479   if(ndigits){
 480     if(verbose){
 481       cout<<"-------------------------------------------------------"<<endl;
 482       cout<<"Number of DIGITS for module "<<mod<<": "<<ndigits<<endl;
 483     }
 484     // Get the coordinates of the module
 485     if(subd==2){
 486       for (Int_t digit=0;digit<ndigits;digit++){
 487         ssdone[digit]=0;
 488         pair[digit]=0;
 489       }
 490     }
 491     for (Int_t digit=0;digit<ndigits;digit++) {
 492       digs = (AliITSdigit*)ITSdigits->UncheckedAt(digit);
 493       Int_t iz=digs->GetCoord1();  // cell number z
 494       Int_t ix=digs->GetCoord2();  // cell number x
 495       // Get local coordinates of the element
 496       if(subd<2){
 497         seg->DetToLocal(ix,iz,lcoor[0],lcoor[2]);
 498       }
 499       else{
 500         // SSD: if iz==0 ---> N side; if iz==1 P side
 501         if(ssdone[digit]==0){
 502           ssdone[digit]=1;
 503           pair[digit]=-1;
 504           Bool_t pside=(iz==1);
 505           Bool_t impaired=kTRUE;
 506           Int_t pstrip=0;
 507           Int_t nstrip=0;
 508           if(pside)pstrip=ix;
 509           if(!pside)nstrip=ix;
 510           for(Int_t digi2=0;digi2<ndigits;digi2++){
 511             if(ssdone[digi2]==0 && impaired){
 512               AliITSdigitSSD *dig2=(AliITSdigitSSD*)ITSdigits->UncheckedAt(digi2);
 513               if(dig2->GetCoord1() != iz && dig2->GetTrack(0)==digs->GetTrack(0) && dig2->GetTrack(0)>0){
 514                 ssdone[digi2]=2;
 515                 pair[digit]=digi2;
 516                 if(pside)nstrip=dig2->GetCoord2();
 517                 if(!pside)pstrip=dig2->GetCoord2();
 518                 impaired=kFALSE;
 519               }
 520             }
 521           }
 522           if(!impaired)seg->GetPadCxz(pstrip,nstrip,lcoor[0],lcoor[2]);
 523         }
 524       }
 525       if(subd<2 || (subd==2 && ssdone[digit]==1)){
 526         Int_t coor1=digs->GetCoord1();
 527         Int_t coor2=digs->GetCoord2();
 528         Int_t tra0=digs->GetTrack(0);
 529         if(verbose){
 530           cout<<"digit # "<<digit<<" fCoord1= "<<coor1<<" fCoord2= ";
 531           cout<<coor2<<" track "<<tra0<<" "<<endl;
 532           if(subd<2)cout<<"local coordinates -- x="<<lcoor[0]<<", z=";
 533           cout<<lcoor[2]<<endl;
 534           if(subd==2){
 535             if(pair[digit]==-1){
 536               cout<<"(digit  not paired)"<<endl;
 537             }
 538             else {
 539               cout<<"local coordinates -- x="<<lcoor[0]<<", z="<<lcoor[2];
 540               cout<<endl;
 541               Int_t dtmp=pair[digit];
 542               AliITSdigitSSD *dig2=(AliITSdigitSSD*)ITSdigits->UncheckedAt(dtmp);
 543               Int_t coor1b=dig2->GetCoord1();
 544               Int_t coor2b=dig2->GetCoord2();
 545               Int_t tra0b=dig2->GetTrack(0);
 546               cout<<"(digit paired with digit #"<<dtmp<<endl;
 547               cout<<"with fCoord1= "<<coor1b<<" fCoord2= "<<coor2b;
 548               cout<<" track "<<tra0b<<")"<<endl;
 549             }
 550           }
 551         }
 552         if(subd<2 || (subd==2 && pair[digit]!=-1)){
 553           // Global coordinates of the element
 554           //SDD and SPD use cm, SSD microns (GetPadCxz)
 555           if(subd==2)for(Int_t j=0;j<3;j++)lcoor[j]=lcoor[j]/10000.;
 556           lcoor[1]=0.;
 557           geom->LtoG(mod,lcoor,gcoor);  // global coord. in cm
 558           ragdig=sqrt(gcoor[0]*gcoor[0]+gcoor[1]*gcoor[1]);
 559           if(verbose){
 560              cout<<"global coordinates "<<gcoor[0]<<" "<<gcoor[1];
 561              cout<<" "<<gcoor[2]<<" Radius "<<ragdig<<endl;
 562           }
 563           //Fill histograms
 564           TH2F *bidi = (TH2F*)histos.At(subd*9);
 565           TH1F *uni = (TH1F*)histos.At(3+subd*9);
 566           bidi->Fill(gcoor[0],gcoor[1]);
 567           uni->Fill(gcoor[2]);
 568         }
 569       }
 570     } // loop on digits for this module
 571   } // if(ndigits>0....
 572 }