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