ITS/ITSgeoplot.C

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