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