Updated version of the Bari code to work with the HEAD. A new test macros has also...
[u/mrichter/AliRoot.git] / ITS / ITSReadPlotData.C
CommitLineData
4913f68a 1Int_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
faafc0eb 120 dtype = (ID[1] - 1) / 2;
4913f68a 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
faafc0eb 162 Text_t msg[250];
163 Float_t xm = 0.0, ym = 0.0, zm = 0.0;
4913f68a 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 }
43c39b72 182 else
183 cout << "Hits scanned..." << endl;
4913f68a 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 }
43c39b72 192 else
193 cout << "Recpoints scanned..." << endl;
4913f68a 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 }
43c39b72 202 else
203 cout << "Digits scanned..." << endl;
4913f68a 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();
43c39b72 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 }
4913f68a 259 }
260
261 cout << "Done. Goodbye" << endl;
262 return;
263}
264
265Int_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
299Int_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
329Int_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
faafc0eb 348 if(dtype==2){
349 Int_t layer, ladder, detec;
350 gm->GetModuleId(ID,layer,ladder,detec);
351 seg->SetLayer(layer);
352 }
353
4913f68a 354 if (!digits_num)
355 return 0;
356 else {
43c39b72 357 cout << "Digits to scan: " << digits_num << endl;
4913f68a 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++) {
4913f68a 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)
43c39b72 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 ***************************************
4913f68a 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]);
43c39b72 406 X[j] /= 10000.0; // convert microns to cm
407 Z[j] /= 10000.0; // convert microns to cm
4913f68a 408 }
409 }
4913f68a 410 }
411 return digits_num;
412}
413