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