]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - ITS/oldmacros/AliITSReadPlotData.C
Fixes for wrong use of const causing PW.CAST_TO_QUALIFIED_TYPE defect in Coverity
[u/mrichter/AliRoot.git] / ITS / oldmacros / AliITSReadPlotData.C
... / ...
CommitLineData
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
30Int_t GetModuleHits (TObject* its, Int_t ID, Float_t*& X, Float_t*& Y, Float_t*& Z, Bool_t*& St);
31Int_t GetModuleRecPoints (TObject *its, Int_t ID, Float_t*& X, Float_t*& Z, Int_t isfastpoints);
32Int_t GetModuleDigits(TObject *its, Int_t ID, Int_t dtype, Float_t*& X, Float_t*& Z);
33
34Int_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
303Int_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
337Int_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
383Int_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