ITS/AliITSGeoPlot.C

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