Updated version of the Bari code to work with the HEAD. A new test macros has also...
[u/mrichter/AliRoot.git] / ITS / ITSReadPlotData.C
1 Int_t ITSReadPlotData(char *filename = "galice.root", Int_t evNum = 0) {
2
3         /*********************************************************************
4          *                                                                   *
5          *  Macro used to read hits, digits and recpoints for a module       *
6          *  It draws 3 TH2Fs where stores the 2-D coordinates of these       *
7          *  data for a specified module (for which the macro asks when       *
8          *  it starts, with some checks to avoid wrong detector coords.      *
9          *                                                                   *
10          *  Only a little 'experimental' warning:                            *
11          *  with all the tests I have made, I wasn't able to plot the        *
12          *  digits fo the 5th layer...                                       *
13          *  I skipped this problem with an alert to beware th user, while    *
14          *  in this situation the macro plots only recpoints and hits        *
15          *                                                                   *
16          *  Author: Alberto Pulvirenti                                       *
17          *                                                                   *
18          *********************************************************************/
19
20         
21         extern Int_t GetModuleHits (TObject* its, Int_t ID[0], Float_t*& X, Float_t*& Y, Float_t*& Z, Bool_t*& St);
22         extern Int_t GetModuleRecPoints (TObject *its, Int_t ID[0], Float_t*& X, Float_t*& Z);
23         extern Int_t GetModuleDigits(TObject *its, Int_t ID[0], Int_t dtype, Float_t*& X, Float_t*& Z);
24         extern void  AssignCoords(TArrayI *ID);
25         
26         // First of all, here are put some variable declarations
27         // that are useful in the following part:
28         Int_t nparticles; // number of particles
29         // ITS module coordinates [layer = 1, ladder = 2, det = 3] and absolute ID[0] of module [0]
30         TArrayI ID(4);
31         Int_t nmodules, dtype; // Total number of modules and module type (SSD, SPD, SDD)
32         Float_t *x = 0, *y = 0, *z = 0; // Arrays where to store read coords
33         Bool_t *St = 0; // Status of the track (hits only)
34
35         // It's necessary to load the gAlice shared libs
36         // if they aren't already stored in memory...
37         if (gClassTable->GetID("AliRun") < 0) {
38         gROOT->LoadMacro("loadlibs.C");
39                 loadlibs();
40         }
41   // Anyway, this macro needs to read a gAlice file, so it
42   // clears the gAlice object if there is already one in memory...
43   else {
44                 if(gAlice){
45                         delete gAlice;
46                         gAlice = 0;
47                 }
48         }
49
50         // Now is opened the Root input file containing Geometry, Kine and Hits
51   // by default its name must be "galice.root".
52   // When the file is opened, its contens are shown.
53         TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(filename);
54         if (!file) file = new TFile(filename);
55                 file->ls();
56         
57         // Then, the macro gets the AliRun object from file.
58         // If this object is not present, an error occurs
59         // and the execution is stopped.
60         // Anyway, this operation needs some time,
61         // don't worry about an apparent absence of output and actions...
62         cout << "\nSearching in '" << filename << "' for an AliRun object ... " << flush;
63         gAlice = (AliRun*)file->Get("gAlice");
64         if (gAlice)
65                 cout << "FOUND!" << endl;
66         else {
67                 cout<<"NOT FOUND! The Macro can't continue!" << endl;
68                 return 0;
69         }
70         
71         // Then, the macro selects the event number specified. Default is 0.
72         nparticles = gAlice->GetEvent(evNum);
73         cout << "\nNumber of particles   = " << nparticles << endl;
74         if (!nparticles) {
75                 cout << "With no particles I can't do much... Goodbye!" << endl;
76                 return 0;
77         }
78         
79         // The next step is filling the ITS from the AliRun object.
80         AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
81   ITS->InitModules(-1, nmodules);
82   cout << "Number of ITS modules = " << nmodules << endl;
83         cout << "\nFilling modules (it takes a while, now)..." << flush;
84         ITS->FillModules(0, 0, nmodules, " ", " ");
85   cout << "DONE!" << endl;
86         AliITSgeom *gm = ITS->GetITSgeom();
87         AliITSDetType *det = ITS->DetType(dtype);       
88         AliITSsegmentation *seg = det->GetSegmentationModel();  
89         
90         for(;;) {
91
92     // Input phase.
93     // The macro asks if the user wants to put a detector ID[0]
94     // or prefers to input layer, ladder and detector.
95     for (Int_t i = 0; i < 4; i++) ID[i] = -1;
96     Int_t answ;
97     do {
98            cout << "\nSelection modes:" << endl;
99                 cout << "1) by layer - ladder - detector numbers" << endl;
100         cout << "2) by unique detector ID" << endl;
101            cout << "0) exit macro" << endl;
102                 cout << "\nEnter your choice: ";
103         cin >> answ;
104     } while (answ < 0 || answ > 2);
105     switch (answ) {
106         case 0:
107                 // input = 0 ---> EXIT
108                 return;
109                 break;
110         case 1:
111                 // input = 1 ---> SELECTION BY COORD
112                 do {
113                                         cout << "\nLayer number [1-6, 0 to exit]: ";
114                                         cin >> ID[1];
115                                         if (!ID[1]) return 0;
116                                 } while (ID[1] < 0 || ID[1] > 6);
117                                 
118                                 // Detector type: 0 --> SPD, 1 --> SDD, 2 --> SSD.
119                                 // Layer 1,2 --> 0 / Layer 3,4 --> 1 / Layer 5,6 --> 2
120                                 dtype = (ID[1] - 1) / 2;
121                                 
122                                 // Once fixed the layer number, the macro calculates the max number
123                                 // for ladder and detector from geometry, and accepts only suitable values.
124                                 do {
125                                         ID[2] = gm->GetNladders(ID[1]);
126                                         cout << "Ladder number [1-" << ID[2] << ", 0 to exit]: ";
127                                         cin >> ID[2];
128                                         if (!ID[2]) return 0;
129                                 } while (ID[2] < 0 || ID[2] > gm->GetNladders(ID[1]));
130                                 do {
131                                         ID[3] = gm->GetNdetectors(ID[1]);
132                                         cout << "Detector number [1-" << ID[3] << ", 0 to exit]: ";
133                                         cin >> ID[3];
134                                         if (!ID[3]) return 0;
135                                 } while (ID[3] < 0 || ID[3] > gm->GetNdetectors(ID[1]));
136                                 break;
137         case 2:
138                 // input = 2 ---> SELECTION BY ID[0]
139                 do {
140                                         ID[0] = gm->GetIndexMax();
141                                         cout << "\n Detector ID number [0-" << ID[0] << ", -1 to exit]: ";
142                                         cin >> ID[0];
143                                         if (ID[0] == -1) return 0;
144                                 } while (ID[0] < 0 || ID[0] > gm->GetIndexMax());
145           break;
146     };
147
148     if (ID[0] == -1)
149                 // If ID[0] = -1 the chioce was by coords, so it's necessary to assign the ID:
150                         ID[0] = gm->GetModuleIndex(ID[1], ID[2], ID[3]);
151                 else {  
152                         // Else we must get the layer, ladder and detector by the ID;
153                         // ...but first we must remember that the ID goest from 0 to NModules - 1!
154                         ID[0]--;
155                         ID[1] = ITS->GetModule(ID[0])->GetLayer();
156                         ID[2] = ITS->GetModule(ID[0])->GetLadder();
157                         ID[3] = ITS->GetModule(ID[0])->GetDet();
158                 }
159                                 
160                 // Defines the histograms inside the `for' cycle, so they are destroyed at the end
161                 // of every read sequqnce, in order to mek another withour segmentation faults
162                 Text_t msg[250];
163                 Float_t xm = 0.0, ym = 0.0, zm = 0.0;
164                 switch (dtype) {
165                         case 0: xm = 1.5; zm = 7.0; break;
166                         case 1: xm = 7.5; zm = 8.0; break;
167                         case 2: xm = 7.5; zm = 4.5; break;
168                 }
169                 sprintf(msg, "Module index=%d lay=%d lad=%d det=%d", ID[0], ID[1], ID[2], ID[3]);
170                 TH2F *hhits = new TH2F("hhits", msg, 500, -xm, xm, 500, -zm, zm);     // Histogram of hits
171                 TH2F *hrecs = new TH2F("hrecs", msg, 500, -xm, xm, 500, -zm, zm);     // Histogram of recpoints
172                 TH2F *hdigits = new TH2F("hdigits", msg, 500, -xm, xm, 500, -zm, zm); // Histogram of digits
173                 
174                 cout << endl;
175                 
176                 // Reads hits...
177                 Int_t hits = GetModuleHits(ITS, ID[0], x, y, z, St);
178                 if (!hits) {
179                         cout << "No hits in module!" << endl;
180                         continue;
181                 }
182                 else
183                         cout << "Hits scanned..." << endl;
184                 for (Int_t i = 0; i < hits; i++) if (!St[i]) hhits->Fill(x[i], z[i]);
185                 
186                 // Reads recpoints...
187                 Int_t recs = GetModuleRecPoints(ITS, ID[0], x, z);
188                 if (!recs) {
189                         cout << "No recpoints in module!" << endl;
190                         continue;
191                 }
192                 else
193                         cout << "Recpoints scanned..." << endl;
194                 for (Int_t i = 0; i < recs; i++) hrecs->Fill(x[i], z[i]);
195                 
196                 // Reads digits...
197                 Int_t digits = GetModuleDigits(ITS, ID[0], dtype, x, z);
198                 if (!digits) {
199                         cout << "No digits in module!" << endl;
200                         //continue;
201                 }
202                 else
203                         cout << "Digits scanned..." << endl;
204                 for (Int_t i = 0; i < digits; i++) hdigits->Fill(x[i], z[i]);
205
206                 // Draws read data...
207                 // HITS -------> red (2) crosses.
208                 // DIGITS -----> green (8) boxes.
209                 // REC-POINTS -> blue (4) St. Andrew's crosses.
210
211                 TCanvas *viewer = new TCanvas("viewer", "Module viewer canvas", 0, 0, 800, 800);
212                 viewer->cd();
213                 
214                 hdigits->SetMarkerStyle(25);
215                 hdigits->SetMarkerColor(8);
216                 hdigits->SetMarkerSize(2);
217                 hdigits->SetStats(kFALSE);
218                 hdigits->SetXTitle("Local X (cm)");
219                 hdigits->SetYTitle("Local Z (cm)");
220                 hdigits->Draw();
221                 
222                 hhits->SetMarkerStyle(5);
223                 hhits->SetMarkerColor(2);
224                 hhits->SetMarkerSize(3);
225                 hhits->SetStats(kFALSE);
226                 hhits->Draw("same");
227         
228                 hrecs->SetMarkerStyle(2);
229                 hrecs->SetMarkerColor(4);
230                 hrecs->SetMarkerSize(3);
231                 hrecs->SetStats(kFALSE);
232                 hrecs->Draw("same");
233                 
234                 TLegend *legend = new TLegend(0.7, 0.8, 0.9, 0.9);
235                 legend->SetMargin(0.2);
236                 legend->AddEntry(hhits, "hits","P");
237                 legend->AddEntry(hrecs, "recpoints","P");
238                 legend->AddEntry(hdigits, "digits","P");
239                 legend->SetTextSizePixels(14);
240                 legend->Draw();
241                 
242                 viewer->Update();
243                 
244                 
245                 Text_t fname[250],ans;
246                 cout << "Do you want to save the current canvas on a file (y/n) ? ";
247                 cin >> ans;
248                 if(ans == 'y' || ans == 'Y') {
249                    cout << "Enter filename: ";
250                    cin >> fname;
251                    TString *control = new TString(fname);
252                    Bool_t ok=control->Contains(".C") || control->Contains(".root") || control->Contains(".ps") || control->Contains(".eps") || control->Contains(".gif");
253                    if(!ok){ 
254                       cout << "File extension is not recognized. The canvas will be saved as Postscript file";
255                       strcat(fname,".ps");
256                    }  
257                    viewer->SaveAs(fname);
258                 }
259         }
260         
261         cout << "Done. Goodbye" << endl;
262         return;
263 }
264
265 Int_t GetModuleHits (TObject* its, Int_t ID, Float_t*& X, Float_t*& Y, Float_t*& Z, Bool_t*& St) {      
266         // First of all, the macro selects the specified module,
267         // then gets the array of hits in it and their number.
268         AliITS *ITS = (AliITS*) its;
269         AliITSmodule *module = ITS->GetModule(ID);
270         TObjArray *hits_array = module->GetHits();
271         Int_t hits_num = hits_array->GetEntriesFast();
272         
273         // Now, if this count returns 0, there's nothing to do,
274         // while, if it doesn't, the first thing to do is dimensioning
275         // the coordinate arrays, and then the loop can start.
276         if (!hits_num)
277                 return 0;
278         else {
279                 if (X) delete [] X;     
280                 if (Y) delete [] Y;
281                 if (Z) delete [] Z;
282                 if (St) delete [] St;
283                 X = new Float_t[hits_num];
284                 Y = new Float_t[hits_num];
285                 Z = new Float_t[hits_num];
286                 St = new Int_t[hits_num];
287         }
288
289         for (Int_t j = 0; j < hits_num; j++) {
290                 AliITShit *hit = (AliITShit*) hits_array->At(j);
291                 X[j]  = hit->GetXL();
292                 Y[j]  = hit->GetYL();
293                 Z[j]  = hit->GetZL();
294                 St[j] = hit->StatusEntering();
295         }
296         return hits_num;
297 }
298
299 Int_t GetModuleRecPoints (TObject *its, Int_t ID, Float_t*& X, Float_t*& Z) {
300         
301         // First of all, the macro selects the specified module,
302         // then gets the array of recpoints in it and their number.
303         AliITS *ITS = (AliITS*) its;
304         TTree *TR = gAlice->TreeR();
305         ITS->ResetRecPoints();
306         TR->GetEvent(ID);
307         TClonesArray *recs_array = ITS->RecPoints();
308         Int_t recs_num = recs_array->GetEntries();
309
310         // Now, if this count returns 0, there's nothing to do,
311         // while, if it doesn't, the first thing to do is dimensioning
312         // the coordinate and energy loss arrays, and then the loop can start.
313         if (!recs_num)
314                 return 0;
315         else {
316                 if (X) delete [] X;     
317                 if (Z) delete [] Z;
318                 X = new Float_t[recs_num];
319                 Z = new Float_t[recs_num];
320         }
321         for(Int_t j = 0; j < recs_num; j++) {
322                 AliITSRecPoint *recp = (AliITSRecPoint*)recs_array->At(j);
323            X[j] = recp->GetX();
324                 Z[j] = recp->GetZ();
325         }
326         return recs_num;        
327 }
328
329 Int_t GetModuleDigits(TObject *its, Int_t ID, Int_t dtype, Float_t*& X, Float_t*& Z) {
330
331         // First of all, the macro selects the specified module,
332         // then gets the array of recpoints in it and their number.
333         AliITS *ITS = (AliITS*)its;
334         TTree *TD = gAlice->TreeD();
335         ITS->ResetDigits();
336         TD->GetEvent(ID);
337         TClonesArray *digits_array = ITS->DigitsAddress(dtype);
338         AliITSgeom *gm = ITS->GetITSgeom();     
339         AliITSDetType *det = ITS->DetType(dtype);       
340         AliITSsegmentation *seg = det->GetSegmentationModel();  
341         TArrayI ssdone(5000);  // used to match p and n side digits of strips
342         TArrayI pair(5000);    // as above      
343         Int_t digits_num = digits_array->GetEntries();
344         // Now, if this count returns 0, there's nothing to do,
345         // while, if it doesn't, the first thing to do is dimensioning
346         // the coordinate and energy loss arrays, and then the loop can start.
347
348         if(dtype==2){
349            Int_t layer, ladder, detec;
350            gm->GetModuleId(ID,layer,ladder,detec);
351            seg->SetLayer(layer);
352         }
353
354         if (!digits_num)
355                 return 0;
356         else {
357                 cout << "Digits to scan: " << digits_num << endl;
358                 if (X) delete [] X;                     
359                 if (Z) delete [] Z;
360                 X = new Float_t[digits_num];            
361                 Z = new Float_t[digits_num];
362         }
363         
364         // Get the coordinates of the module
365         if (dtype == 2) {
366                 for (Int_t j = 0; j < digits_num; j++) {
367                         ssdone[j] = 0;                  
368                         pair[j] = 0;
369                 }
370         }
371   for (Int_t j = 0; j < digits_num; j++) {
372                 digit = (AliITSdigit*)digits_array->UncheckedAt(j);
373                 Int_t iz=digit->fCoord1;  // cell number z
374                 Int_t ix=digit->fCoord2;  // cell number x
375     // Get local coordinates of the element (microns)
376                 // ******************************* PARTE CORRETTA ***************************************
377                 if(dtype < 2) {
378                         Float_t xx, zz; // aggiunta
379         seg->DetToLocal(ix, iz, xx, zz);
380                         X[j] = xx; // aggiunta
381                         Z[j] = zz; // aggiunta
382                 }
383                 // ******************************* FINE PARTE CORRETTA ***************************************
384     else {
385                         // SSD: if iz==0 ---> N side; if iz==1 P side
386       if (ssdone[j] == 0) {
387                                 ssdone[j]=1;
388                                 pair[j]=-1;
389                                 Bool_t pside = (iz == 1);
390                                 Bool_t impaired = kTRUE;
391                                 Int_t pstrip = 0;
392                                 Int_t nstrip = 0;
393                                 if (pside) pstrip = ix; else nstrip = ix;
394                                 for (Int_t k = 0; k < digits_num; k++) {
395                                         if (ssdone[k] == 0 && impaired) {
396                                                 AliITSdigitSSD *sdigit=(AliITSdigitSSD*)digits_array->UncheckedAt(k);
397                                                 if (sdigit->fCoord1 != iz && sdigit->GetTracks()[0] == digit->GetTracks()[0]) {
398                                                         ssdone[k]=2;
399                                                         pair[j]=k;
400                                                         if (pside) nstrip = sdigit->fCoord2; else pstrip = sdigit->fCoord2;
401                                                         impaired=kFALSE;
402                                                 }
403                                         }
404                                 }
405         if (!impaired) seg->GetPadCxz(pstrip, nstrip, X[j], Z[j]);
406                                 X[j] /= 10000.0;  // convert microns to cm
407                                 Z[j] /= 10000.0;  // convert microns to cm
408                         }
409                 }
410         }
411         return digits_num;
412
413