bc3a469a419046dc19179182fd01fb68e0bf922e
[u/mrichter/AliRoot.git] / ITS / ITSgeoplot.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 ITSgeoplot (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->GetTracks()[0]==digs->GetTracks()[0] && dig2->GetTracks()[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->GetTracks()[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->GetTracks()[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