Vector renamed to avoid re-definition
[u/mrichter/AliRoot.git] / ITS / ITSOccupancy.C
CommitLineData
2bba24f4 1/******************************************************************************
2
3 "ITSOccupancy.C"
4
5 this macro calculates the mean occupancy of each ITS layer,
6 making also a distribution in function of the z-value of the
7 "fired" digits for each layer
8
9 Because of evaluation problems, it is necessary to permit the
10 analysis for both hits, digits and recpoints, so it's provided an
11 option data which must contain
12 - "H" for hits
13 - "D" for digits
14 - "R" for recpoints
15 but it's possible to use several option in the same time
16 (like "HD", "HRD", and so on...)
17
18 to avoid the problem of drawing unrequested windows, the graph are
19 plotted only if it is explictly requested with the option "P", that
20 can be inserted together to the ones described above.
21
22 It is also possible to compile this macro, to execute it faster.
23 To do this, it is necessary to:
24 1) make sure that de "#define __COMPILED__" instruction below is UNcommented
25 and the "#define __NOCOMPILED__" instruction is commented (remember that
26 if you comment or uncomment both these instructions, the macro can't work)
27 2) type, at the root prompt, the instruction
28 gSystem->SetIncludePath("-I- -I$ALICE_ROOT/include -I$ALICE_ROOT/ITS -g")
29 to adjoin the ALICE include libraries to the main include directory
30 3) load the function via the ".L ITSOccupancy.C++" statement (only the first
31 time, because the first time the macro must be compiled). Frome the
32 second time you use the macro, you must only load it via the
33 ".L ITSOccupancy.C+" instruction (note the number of '+' signs in each case
34 4) call the function with only its name, es: ITSOccupancy("HP", "galice.root")
35
36 NOTE 1: option interpretation is case-insensitive ("h" eqv "H", etc.)
37 NOTE 2: if you insert a filename as argument, it must be inserted only the name
38 (the macro attaches the extension ".root" to it, every time)
39
40 Author: Alberto Pulvirenti
41******************************************************************************/
42
43#define __COMPILED__
44//#define __NOCOMPILED__
45
46#ifdef __COMPILED__
47 #include <iostream.h>
48 #include <TROOT.h>
49 #include <TArrayI.h>
50 #include <TCanvas.h>
51 #include <TClassTable.h>
52 #include <TClonesArray.h>
53 #include <TFile.h>
54 #include <TObject.h>
55 #include <TObjArray.h>
56 #include <TTree.h>
57 #include <TMath.h>
58 #include <TString.h>
59 #include <TH1.h>
60 #include <AliRun.h>
61 #include <AliITS.h>
62 #include <AliITSgeom.h>
63 #include <AliITSDetType.h>
64 #include <AliITSRecPoint.h>
65 #include <AliITSdigit.h>
66 #include <AliITShit.h>
67 #include <AliITSmodule.h>
68 #include <AliITSsegmentation.h>
69 #include <AliITSsegmentationSPD.h>
70 #include <AliITSsegmentationSDD.h>
71 #include <AliITSsegmentationSSD.h>
72#endif
73
74Int_t ITSOccupancy(char* opt, char *name = "galice", Int_t evNum = 0) {
75
76 extern void GetModuleDigits(TObject *its, Int_t ID, Int_t dtype, TH1D *counter);
77 extern void GetModuleHits (TObject* its, Int_t ID, TH1D *counter);
78 extern void GetModuleRecPoints (TObject *its, Int_t ID, TH1D *counter);
79 extern Int_t DType(Int_t layer);
80
81 // Evaluating options
82 TString option(opt);
83 Bool_t HITS = (option.Contains("h") || option.Contains("H"));
84 Bool_t DIGS = (option.Contains("d") || option.Contains("D"));
85 Bool_t RECS = (option.Contains("r") || option.Contains("R"));
86 Bool_t PLOT = (option.Contains("p") || option.Contains("P"));
87
88 Text_t filename[100];
89 char gifH[100], gifD[100], gifR[100];
90 strcpy (filename, name);
91 strcpy (gifH, name);
92 strcpy (gifD, name);
93 strcpy (gifR, name);
94 strcat (filename, ".root");
95 strcat (gifH, "H.gif");
96 strcat (gifD, "D.gif");
97 strcat (gifR, "R.gif");
98
99 if (!HITS && !DIGS && !RECS) {
100 cout << "\n\n--- NO DATA-TYPE SPECIFIED!!! ---\n\n";
101 return 0;
102 }
103
104 // --------------------------------
105 // 1. ALICE objects setting phase
106 // --------------------------------
107
108 // Load the gAlice shared libs if not already in memory
109#ifdef __NOCOMPILED__
110 if (gClassTable->GetID("AliRun") < 0) {
111 gROOT->LoadMacro("loadlibs.C");
112 loadlibs();
113 }
114#endif
115 if(gAlice){
116 delete gAlice;
117 gAlice=0;
118 }
119 // Open the Root input file containing the AliRun data
120 // (default name is "galice.root")
121 TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(filename);
122 if (!file) file = new TFile(filename);
123 file->ls();
124 // Get the AliRun object from file (if there's none, the macro is stopped)
125 cout << "\nSearching in '" << filename << "'" << endl;
126 gAlice = (AliRun*)file->Get("gAlice");
127 if (gAlice)
128 cout << "Ok, I found an AliRun object..." << endl;
129 else {
130 cout<<"Sorry, there isn't any AliRun object here..." << endl;
131 return 0;
132 }
133 // Select the event number specified. Default is 0.
134 Int_t nparticles = gAlice->GetEvent(evNum);
135 cout << "\nNumber of particles = " << nparticles << endl;
136 if (!nparticles) {
137 cout << "With no particles I can't do much..." << endl;
138 return 0;
139 }
140 // Initialize the ITS and its geometry
141 AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
142 AliITSgeom *gm = ITS->GetITSgeom();
143 // Fill the AliITSmodule objects (only for hits)
144 Int_t nmodules;
145 if (HITS) {
146 ITS->InitModules(-1, nmodules);
147 cout << "Number of ITS modules = " << nmodules << endl;
148 cout << "\nFilling modules (it takes a while, now)..." << flush;
149 ITS->FillModules(0, 0, nmodules, " ", " ");
150 cout << "DONE!" << endl;
151 }
152
153 // -------------------------------------------------
154 // B. Variables and objects definition and setting
155 // -------------------------------------------------
156
157 Double_t zmax[6] = {16.5, 16.5, 22.2, 29.7, 45.1, 50.8}; // max z measured (cm)
158 Double_t binw[] = {2., 2., 2., 2., 5., 5.}; // bin width for TH1Ds (cm)
159
160 Double_t tot[6]; // total number of channels
161
162 // Histos for plotting hits [1], digits [2] and points [3], and divisor [0]
163 TH1D *hist[4][6];
164
165 // Histograms initialization:
166 // Using the AliITSsegmentation object, here is deduced
167 // the maximum z that can be measurable for each layer, simply
168 // obtaining the globar coordinates of the first digit of the first module
169 // of each layer (that is placed at the maximum z).
170 for (Int_t i = 0; i < 6; i++) {
171 AliITSDetType *det = ITS->DetType(DType(i));
172 AliITSsegmentation *seg = det->GetSegmentationModel();
173 tot[i] = gm->GetNladders(i+1) * gm->GetNdetectors(i+1) * seg->GetNPads();
174 cout.setf(ios::fixed);
175 cout << tot[i] << endl;
23a8029f 176 Text_t hname[7];
2bba24f4 177 for (Int_t j = 0; j < 4; j++) {
23a8029f 178 hname[0] = 'h';
179 hname[1] = 'i';
180 hname[2] = 's';
181 hname[3] = 't';
182 hname[4] = '0' + j;
183 hname[5] = '0' + i;
184 hname[6] = '\0';
2bba24f4 185 Int_t nbins = (Int_t)(2. * zmax[i] / binw[i]);
23a8029f 186 hist[j][i] = new TH1D(hname, "histo", nbins, -zmax[i], zmax[i]);
2bba24f4 187 }
188 }
189 // Here is put the value of the density for each bin of each layer's TH1D
190 for (Int_t i = 0; i < 6; i++) {
191 for (Int_t j = 0; j < hist[0][i]->GetNbinsX(); j++) {
192 Double_t z = hist[0][i]->GetBinLowEdge(j);
193 Double_t bincontent = tot[i] / (zmax[i] * 2.0) * hist[0][i]->GetBinWidth(j);
194 hist[0][i]->Fill(z + 0.1, bincontent / 100.0);
195 z += hist[0][i]->GetBinWidth(j);
196 }
197 }
198
199 // --------------------------------------------------------------
200 // Counting the hits, digits and recpoints contained in the ITS
201 // --------------------------------------------------------------
202
203 for (Int_t L = 0; L < 6; L++) {
204
205 // To avoid two nested loops, are calculated
206 // the ID of the first and last module of the L
207 // (remember the L goest from 0 to 5 (not from 1 to 6)
208 Int_t first, last;
209 first = gm->GetModuleIndex(L + 1, 1, 1);
210 last = gm->GetModuleIndex(L + 1, gm->GetNladders(L + 1), gm->GetNdetectors(L + 1));
211
212 // Start loop on modules
213 cout << "Examinating layer " << L + 1 << " ... " << flush;
214 for (Int_t ID = first; ID <= last; ID++) {
215
216 // Hits analysis (if requested)
217 if (HITS)
218 GetModuleHits(ITS, ID, hist[2][L]);
219
220 // Digits analysis (if requested)
221 if (DIGS)
222 GetModuleDigits(ITS, ID, DType(L), hist[1][L]);
223
224 // Recpoints analysis (if requested)
225 if (RECS)
226 GetModuleRecPoints(ITS, ID, hist[3][L]);
227 }
228
229 hist[1][L]->Divide(hist[0][L]);
230 hist[2][L]->Divide(hist[0][L]);
231 hist[3][L]->Divide(hist[0][L]);
232 cout << "DONE" << endl;
233 }
234
235 // Write on the screen the total means
236 cout << endl;
237 cout << "********* MEAN OCCUPANCIES *********" << endl;
238 for (Int_t L = 0; L < 6; L++) {
239 Double_t mean[3];
240 cout << " LAYER " << L << ": " << endl;
241 for (Int_t i = 0; i < 3; i++) {
242 mean[i] = hist[i+1][L]->GetSumOfWeights() / hist[i+1][L]->GetNbinsX();
243 Text_t title[50];
244 sprintf(title, " - Layer %d, mean %4.2f - ", L + 1, mean[i]);
245 hist[i+1][L]->SetTitle(title);
246 cout.setf(ios::fixed);
247 cout.precision(3);
248 if (HITS && i == 1) {
249 cout << " hits occupancy = ";
250 if (mean[i] < 10.0) cout << ' ';
251 cout << mean[i] << "%" << endl;
252 }
253 else if (DIGS && i == 0) {
254 cout << " digits occupancy = ";
255 if (mean[i] < 10.0) cout << ' ';
256 cout << mean[i] << "%" << endl;
257 }
258 if (RECS && i == 2) {
259 cout << " recpts occupancy = ";
260 if (mean[i] < 10.0) cout << ' ';
261 cout << mean[i] << "%" << endl;
262 }
263 }
264 cout << "------------------------------------" << endl;
265 }
266
267 // ----------------------
268 // Plots (if requested)
269 // ----------------------
270 if (!PLOT) return 1;
271
272 TCanvas *view[4];
273 if (HITS) {
274 view[2] = new TCanvas("viewH", "Occupancy view (HITS)", 0, 0, 1050, 700);
275 view[2]->Divide(3,2, 0.001, 0.001);
276 }
277 if (DIGS) {
278 view[1] = new TCanvas("viewD", "Occupancy view (DIGITS)", 20, 40, 1050, 700);
279 view[1]->Divide(3,2, 0.001, 0.001);
280 }
281 if (RECS) {
282 view[3] = new TCanvas("viewR", "Occupancy view (RECPOINTS)", 40, 60, 1050, 700);
283 view[3]->Divide(3,2, 0.001, 0.001);
284 }
285
286 for (Int_t L = 0; L < 6; L++) {
287 for (Int_t i = 1; i <= 3; i++) {
288 if (DIGS && i == 1)
289 hist[i][L]->SetFillColor(kBlue);
290 else if (HITS && i == 2)
291 hist[i][L]->SetFillColor(kRed);
292 else if (RECS && i == 3)
293 hist[i][L]->SetFillColor(kGreen);
294
295 if ((HITS && i == 2) || (DIGS && i == 1) || (RECS && i == 3)) {
296 view[i]->cd(L+1);
297 hist[i][L]->SetStats(kFALSE);
298 hist[i][L]->Draw();
299 hist[i][L]->GetXaxis()->SetTitle("zeta");
300 hist[i][L]->GetYaxis()->SetTitle("%");
301 view[i]->Update();
302 }
303 }
304 }
305
306 if (HITS) view[2]->SaveAs(gifH);
307 if (DIGS) view[1]->SaveAs(gifD);
308 if (RECS) view[3]->SaveAs(gifR);
309
310 return 1;
311}
312
313void GetModuleDigits(TObject *its, Int_t ID, Int_t dtype, TH1D *counter) {
314 // First of all, the macro selects the specified module,
315 // then gets the array of recpoints in it and their number.
316 AliITS *ITS = (AliITS*)its;
317 TTree *TD = gAlice->TreeD();
318 ITS->ResetDigits();
319 TD->GetEvent(ID);
320 TClonesArray *digits_array = ITS->DigitsAddress(dtype);
321 AliITSgeom *gm = ITS->GetITSgeom();
322 AliITSDetType *det = ITS->DetType(dtype);
323 AliITSsegmentation *seg = det->GetSegmentationModel();
324 Int_t digits_num = digits_array->GetEntries();
325
326 // Get the coordinates of the module
327 for (Int_t j = 0; j < digits_num; j++) {
328 Float_t lx[] = {0.0,0.0,0.0}, gx[] = {0.0,0.0,0.0};
329 AliITSdigit *digit = (AliITSdigit*)digits_array->UncheckedAt(j);
330 Int_t iz=digit->fCoord1; // cell number z
331 Int_t ix=digit->fCoord2; // cell number x
332 // Get local coordinates of the element (microns)
333 seg->GetPadCxz(ix, iz, lx[0], lx[2]);
334 if (dtype != 1) {
335 // !!!THIS CONVERSION TO HIT LRS SHOULD BE REMOVED AS SOON AS THE CODE IS FIXED
336 if (!dtype) {
337 lx[0] -= seg->Dx() / 2.0;
338 lx[2] -= seg->Dz() / 2.0;
339 }
340 lx[0] /= 10000.0; // convert from microns to cm (if not SDD)
341 lx[2] /= 10000.0; // convert from microns to cm (if not SDD)
342 }
343 gm->LtoG(ID, lx, gx);
344 counter->Fill(gx[2]);
345 }
346}
347
348/*void GetModuleDigits(TObject *its, Int_t ID, Int_t dtype, TH1D *counter) {
349 // First of all, the macro selects the specified module,
350 // then gets the array of recpoints in it and their number.
351 AliITS *ITS = (AliITS*)its;
352 TTree *TD = gAlice->TreeD();
353 ITS->ResetDigits();
354 TD->GetEvent(ID);
355 TClonesArray *digits_array = ITS->DigitsAddress(dtype);
356 AliITSgeom *gm = ITS->GetITSgeom();
357 AliITSDetType *det = ITS->DetType(dtype);
358 AliITSsegmentation *seg = det->GetSegmentationModel();
359 TArrayI ssdone(5000); // used to match p and n side digits of strips
360 TArrayI pair(5000); // as above
361 Int_t digits_num = digits_array->GetEntries();
362
363 // Get the coordinates of the module
364 if (dtype == 2) {
365 for (Int_t j = 0; j < digits_num; j++) {
366 ssdone[j] = 0;
367 pair[j] = 0;
368 }
369 }
370 for (Int_t j = 0; j < digits_num; j++) {
371 Double_t lx[] = {0.0,0.0,0.0}, gx[] = {0.0,0.0,0.0};
372 AliITSdigit *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 if(dtype < 2) {
377 seg->GetPadCxz(ix, iz, lx[0], lx[2]);
378 if (dtype == 0) {
379 // !!!THIS CONVERSION TO HIT LRS SHOULD BE REMOVED AS SOON AS THE CODE IS FIXED
380 lx[0] -= seg->Dx() / 2.0;
381 lx[2] -= seg->Dz() / 2.0;
382 lx[0] /= 10000.0; // convert from microns to cm (SPD)
383 lx[2] /= 10000.0; // convert from microns to cm (SPD)
384 }
385 gm->LtoG(ID, lx, gx);
386 counter->Fill(gx[2]);
387 }
388 else {
389 // SSD: if iz==0 ---> N side; if iz==1 P side
390 if (ssdone[j] == 0) {
391 ssdone[j]=1;
392 pair[j]=-1;
393 Bool_t pside = (iz == 1);
394 Bool_t impaired = kTRUE;
395 Int_t pstrip = 0;
396 Int_t nstrip = 0;
397 if (pside) pstrip = ix; else nstrip = ix;
398 for (Int_t k = 0; k < digits_num; k++) {
399 if (ssdone[k] == 0 && impaired) {
400 AliITSdigitSSD *sdigit=(AliITSdigitSSD*)digits_array->UncheckedAt(k);
401 if (sdigit->fCoord1 != iz && sdigit->GetTracks()[0] == digit->GetTracks()[0]) {
402 ssdone[k]=2;
403 pair[j]=k;
404 if (pside) nstrip = sdigit->fCoord2; else pstrip = sdigit->fCoord2;
405 impaired=kFALSE;
406 }
407 }
408 }
409 if (!impaired) {
410 seg->GetPadCxz(pstrip, nstrip, lx[0], lx[2]);
411 lx[0] /= 10000.0;
412 lx[2] /= 10000.0;
413 gm->LtoG(ID, lx, gx);
414 counter->Fill(gx[2]);
415 }
416 }
417 }
418 }
419}*/
420
421void GetModuleHits (TObject* its, Int_t ID, TH1D *counter) {
422 // First of all, the macro selects the specified module,
423 // then gets the array of hits in it and their number.
424 AliITS *ITS = (AliITS*) its;
425 AliITSmodule *module = ITS->GetModule(ID);
426 TObjArray *hits_array = module->GetHits();
427 Int_t hits_num = hits_array->GetEntriesFast();
428 // next, fills the counters with the hits coordinates' calcula
429 for (Int_t j = 0; j < hits_num; j++) {
430 AliITShit *hit = (AliITShit*) hits_array->At(j);
431 Double_t Z = hit->GetZG();
432 if (!hit->StatusEntering()) counter->Fill(Z);
433 }
434}
435
436void GetModuleRecPoints (TObject *its, Int_t ID, TH1D *counter) {
437 // First of all, the macro selects the specified module,
438 // then gets the array of recpoints in it and their number.
439 AliITS *ITS = (AliITS*) its;
440 AliITSgeom *gm = ITS->GetITSgeom();
441 TTree *TR = gAlice->TreeR();
442 if (!TR || TR->GetEntries() == 0) {
443 cout << "The ITS wasn't clustered..." << endl;
444 return;
445 }
446 ITS->ResetRecPoints();
447 TR->GetEvent(ID);
448 TClonesArray *recs_array = ITS->RecPoints();
449 Int_t recs_num = recs_array->GetEntries();
450
451 for(Int_t j = 0; j < recs_num; j++) {
452 Double_t lx[3] = {0.,0.,0.}, gx[3] = {0.,0.,0.};
453 AliITSRecPoint *recp = (AliITSRecPoint*)recs_array->At(j);
454 lx[0] = recp->GetX();
455 lx[1] = 0.0;
456 lx[2] = recp->GetZ();
457 gm->LtoG(ID, lx, gx);
458 counter->Fill(gx[2]);
459 }
460}
461
462Int_t DType(Int_t layer) {
463 if (layer == 0 || layer == 1)
464 return 0;
465 else if (layer == 2 || layer == 3)
466 return 1;
467 else
468 return 2;
469}