]>
Commit | Line | Data |
---|---|---|
9f3c053c | 1 | /** Welcome! |
2 | * | |
3 | * This macro is intended for the following tasks: | |
4 | * 1. find bad (switched off/dead/noisy/strange) cells; | |
5 | * 2. find out the status/quality of analysed data (missing RCUs, run quality); | |
6 | * 3. find the extent of problems related to bad cells: required for | |
7 | * systematics estimation for a physical quantity related to a cluster. | |
8 | * | |
9 | * For impatient, change the "infile" line below and execute | |
10 | * root -l getCellsRunsQA.C | |
11 | * | |
12 | * For curios, continue reading getCellsRunsQA() code: it is self-documented. | |
13 | * | |
14 | * The macro options are tuned for a user (and pp runs), and in most cases no user | |
15 | * intervension is necessary. Still, it is likely that you will have to edit | |
16 | * nexc/excells[] in the parameters section below and run the macro several times. | |
17 | * Consult with the output from this macro. | |
18 | * In case of PbPb runs, a small modification is necessary: | |
19 | * 1) change ExcludeSmallRuns() line by putting a smaller number, e.g. 5-10k events; | |
20 | * 2) change FindDeadNoisyCellsPerRun() factor thresholds to a more narrow region, e.g. 0.07-2; | |
21 | * 3) probably, change fit region in FitPi0(). | |
22 | * Also, do not forget to adjust cluster cut for pi0s in AliCaloCellsQA. The value 2.5GeV | |
23 | * is currently reasonable. | |
24 | * | |
25 | * Generally, a QA expert uncomments all the functions which return (print to stdout) | |
26 | * bad cell candidates and checks them by hand. | |
27 | * | |
28 | * Detector-specific parts require to run this macro with aliroot instead of root, | |
29 | * they are commented in getCellsRunsQA() by default. | |
30 | * | |
31 | * This macro is written as a number of small functions and is designed both | |
32 | * for EMCAL and for PHOS (and DCAL in future). | |
33 | * Drawing options are chosen with respect to PPRStyle() drawing defaults. | |
34 | * | |
35 | * Input: AliCaloCellsQA analysis output. | |
36 | * | |
37 | * TODO: cells time spectra: currently it is not put in use. Seems that time shape fitting is | |
38 | * not trivial due to the presence of parasite peaks (one needs to remove them first?), | |
39 | * and this is a separate issue to think about... | |
40 | * | |
41 | * TODO: some PHOS-specific parts | |
42 | * | |
43 | * Author: Olga Driga (SUBATECH) | |
44 | */ | |
45 | ||
46 | ||
47 | // ========================== Parameters section ============================== | |
48 | ||
49 | // input | |
35ac8df7 | 50 | // char *infile = "LHC11e_cpass1_CellQA_AnyInt.root"; |
51 | char *infile = "LHC11e_cpass1_CellQA_PHI7.root"; | |
52 | ||
53 | Bool_t trendPlotEvents = kFALSE; | |
9f3c053c | 54 | |
55 | // supermodule colors | |
56 | Color_t colorsSM[] = {0,2,3,4}; | |
57 | ||
58 | // cells to exclude from averages calculations and their number: bad cells can | |
59 | // mess up the results, it is suggested to list (known and newly found with | |
60 | // this code) problematic cell candidates explicitly | |
61 | ||
35ac8df7 | 62 | Int_t excells[] = {1}; |
63 | Int_t nexc = 0; | |
9f3c053c | 64 | |
65 | // ====================== End of parameters section =========================== | |
66 | ||
67 | static TFile *gInputFile; | |
68 | ||
69 | void getCellsRunsQAPHOS() | |
70 | { | |
71 | // Entry point for the analysis. | |
72 | ||
73 | gRandom->SetSeed(0); | |
74 | gStyle->SetOptStat(0); | |
75 | gStyle->SetOptFit(111); | |
76 | gStyle->SetPalette(1); | |
77 | gStyle->SetFillColor(10); | |
78 | ||
79 | Printf("Input: %s", infile); | |
80 | gInputFile = TFile::Open(infile); | |
81 | ||
82 | // You may wish to extract and draw a particular histogram from infile. | |
83 | // Here are the examples: | |
84 | // | |
85 | // TH1* hNEventsProcessedPerRun = (TH1*) gInputFile->Get("hNEventsProcessedPerRun"); | |
86 | // hNEventsProcessedPerRun->Draw(); | |
87 | // return; | |
88 | // | |
89 | // TH1* hNTimesInClusterElow = (TH1*) gInputFile->Get("run146686_hCellLocMaxNTimesInClusterElow"); | |
90 | // hNTimesInClusterElow->Draw(); | |
91 | // return; | |
92 | // | |
93 | // TH1* hETotalClusterEhigh = (TH1*) gInputFile->Get("run146686_hCellLocMaxETotalClusterEhigh"); | |
94 | // hETotalClusterEhigh->Draw(); | |
95 | // return; | |
96 | ||
97 | // Draw a random cell spectrum; | |
98 | // 0.2-1 GeV -- fit region; | |
99 | // 3 GeV -- the spectrum is drawn up to this value, -1 = no limit; | |
100 | // last argument -- histogram name to process, possible values are: | |
101 | // hCellAmplitude, hCellAmplitudeEHigh, hCellAmplitudeNonLocMax or fhCellAmplitudeEhighNonLocMax. | |
102 | DrawCell(1+gRandom->Integer(10752), 0.25, 1., 4., "hCellAmplitude"); // common cell region for EMCAL2010 and for PHOS | |
103 | ||
104 | // Draw a random cell time spectrum | |
105 | // DrawCellTime(1+gRandom->Integer(10752)); | |
106 | ||
107 | ||
108 | /* RUN NUMBER SELECTION SECTION | |
109 | * | |
110 | * NOTE: at any point below runs are sorted in chronological order. | |
111 | */ | |
112 | ||
113 | // array with run numbers and their number | |
114 | Int_t runNumbers[10000]; | |
115 | Int_t nruns = 0; | |
116 | ||
117 | // First, fill run numbers ... | |
118 | GetRunNumbers(nruns, runNumbers); | |
119 | Printf("Total number of runs: %i", nruns); | |
120 | ||
121 | // ... draw events distribution ... | |
122 | // (the last argument is number of bins in this distribution) | |
123 | DrawRunsDistribution(nruns, runNumbers, 100); | |
124 | ||
125 | // ... and exclude runs with number of events < 1k. | |
c9365469 | 126 | ExcludeSmallRuns(nruns, runNumbers, 100); |
9f3c053c | 127 | |
128 | // You may wish to exclude particular runs: | |
129 | // Int_t runs2Exclude[] = {111222,333444,555666}; | |
130 | // ExcludeRunNumbers(nruns, runNumbers, 3, runs2Exclude); | |
131 | ||
132 | Printf("Number of runs to be analysed: %i", nruns); | |
133 | ||
134 | // Finally, print a nice table with run index / run number / number of events. | |
135 | PrintRunNumbers(nruns, runNumbers); | |
136 | ||
137 | ||
138 | /* PER RUN BAD CELLS SEARCHING CRITERIA | |
139 | * | |
140 | * Four primary criteria on a per run basis: | |
141 | * 1 and 2: number of times cell was a local maximum in a cluster at low/high energies; | |
142 | * 3 and 4: total cluster energy for a local maximum cell at low/high energies. | |
143 | */ | |
144 | ||
145 | // Extract the histograms with dead/noisy cell candidates per run. | |
146 | // For each of the four criteria, for each run and for each cell: | |
147 | // 1) calculate cell factor = [cell value]/[average over cells]; | |
148 | // 2) mark cell as dead (candidate) if factor <= factorDead (3rd argument below); | |
149 | // 3) mark cell as noisy (candidate) if factor >= factorNoisy (4th argument below). | |
150 | // | |
151 | // Factor thresholds are quite wide by default: | |
152 | // low energy criteria are not very sensitive to them, | |
153 | // while high energy criteria are very sensitive due to limited statistics. | |
154 | // | |
155 | // The function below also draws histograms with factor distributions in all the runs | |
156 | // for all the cells. It may help to take the decision about dead/noisy factor thresholds. | |
157 | // The last two arguments -- number of bins and maximum X axis value for such histograms. | |
158 | // | |
159 | // Output: hBadCellMap[4] -- bad cell candidates per run for each of the four criteria; | |
160 | // X axis -- cell absId, Y axis -- run index, content: 0=not bad, 1=noisy, -1=dead. | |
161 | // | |
162 | ||
c9365469 | 163 | TH2** hBadCellMap = FindDeadNoisyCellsPerRun(nruns, runNumbers, 0.05, 2.5, 200, 10.); |
9f3c053c | 164 | |
165 | // Look for noisy channels only: | |
35ac8df7 | 166 | // TH2** hBadCellMap = FindDeadNoisyCellsPerRun(nruns, runNumbers, -1.0, 2.5, 200, 20.); |
9f3c053c | 167 | |
168 | // Look for dead channels only: | |
169 | // TH2** hBadCellMap = FindDeadNoisyCellsPerRun(nruns, runNumbers, 0.05, 1.e9, 200, 10.); | |
170 | ||
171 | // Print cell numbers suggested for exclusion from averages calculation; | |
172 | // see excells[] array in the parameters section in the beginning of the macro. | |
173 | // It is highly suggested to run this macro several times and add the output | |
174 | // of this function to the parameters section. | |
175 | PrintCellsToExcludeFromAverages(hBadCellMap); | |
176 | ||
177 | // The results can be visualized (second argument -- canvas name): | |
178 | // either per run ... | |
179 | DrawDeadNoisyCellMap(hBadCellMap[0], "hMapRuns_NTimesInClusterElow"); | |
180 | DrawDeadNoisyCellMap(hBadCellMap[1], "hMapRuns_NTimesInClusterEhigh"); | |
181 | // DrawDeadNoisyCellMap(hBadCellMap[2], "hMapRuns_ETotalClusterElow"); | |
182 | // DrawDeadNoisyCellMap(hBadCellMap[3], "hMapRuns_ETotalClusterEhigh"); | |
183 | ||
184 | // ... or per number of events ... | |
185 | // DrawDeadNoisyCellMap((TH2*)ConvertMapRuns2Events(nruns,runNumbers,hBadCellMap[0]), "hMapEvents_NTimesInClusterElow"); | |
186 | // DrawDeadNoisyCellMap((TH2*)ConvertMapRuns2Events(nruns,runNumbers,hBadCellMap[1]), "hMapEvents_NTimesInClusterEhigh"); | |
187 | // DrawDeadNoisyCellMap((TH2*)ConvertMapRuns2Events(nruns,runNumbers,hBadCellMap[2]), "hMapEvents_ETotalClusterElow"); | |
188 | // DrawDeadNoisyCellMap((TH2*)ConvertMapRuns2Events(nruns,runNumbers,hBadCellMap[3]), "hMapEvents_ETotalClusterEhigh"); | |
189 | ||
190 | // ... or we can also draw the percent for each cell being dead/noisy | |
191 | // DrawDeadNoisyCellPercent(nruns, runNumbers, hBadCellMap[0], "hPercent_NTimesInClusterElow"); | |
192 | // DrawDeadNoisyCellPercent(nruns, runNumbers, hBadCellMap[1], "hPercent_NTimesInClusterEhigh"); | |
193 | // DrawDeadNoisyCellPercent(nruns, runNumbers, hBadCellMap[2], "hPercent_ETotalClusterElow"); | |
194 | // DrawDeadNoisyCellPercent(nruns, runNumbers, hBadCellMap[3], "hPercent_ETotalClusterEhigh"); | |
195 | ||
196 | // Our main criteria for analytical finding of bad cells is based on the following. | |
197 | // Note that factor distributions for NTimesInCluster and ETotalCluster are very | |
198 | // similar, both at low and at high energy. Note also that we can say the same | |
199 | // for dead/noisy cell visualizations above: they are very similar. This suggests | |
200 | // that NTimesInCluster and ETotalCluster give the same results. The criteria | |
201 | // NTimesInCluster, however, is more calibration-independent (though energy thresholds | |
202 | // are still calibration-dependent) and thus is more reliable and clear. Thus, we | |
203 | // limit ourselves to NTimesInClusterElow/NTimesInClusterEhigh criteria. | |
204 | // Now, you could probably noted that NTimesInClusterEhigh give more dead | |
205 | // cells than that at low energy. This is an expected statistical effect: we have | |
206 | // much smaller number of clusters at high energy. Consequently, we will not use dead | |
207 | // cell candidates at high energy. | |
208 | // | |
209 | // Finally, we combine candidates from low/high energies and produce one TH2 | |
210 | // histogram which is the primary source of our analytical results. | |
211 | // | |
212 | TH2* hBadCellMapPrimary = CombineDeadNoisyElowEhigh(hBadCellMap[0], hBadCellMap[1]); | |
213 | ||
214 | // Note for PHOS: if you are not happy with NTimesInClusterEhigh results (due to a lack of statistics) | |
215 | // uncomment this line: | |
216 | // hBadCellMapPrimary = hBadCellMap[0]; | |
217 | ||
218 | // Draw everything for combined | |
219 | DrawDeadNoisyCellMap(hBadCellMapPrimary, "Primary_hMapRuns"); | |
35ac8df7 | 220 | // DrawDeadNoisyCellMap((TH2*)ConvertMapRuns2Events(nruns,runNumbers,hBadCellMapPrimary), "Primary_hMapEvents"); |
221 | // DrawDeadNoisyCellPercent(nruns, runNumbers, hBadCellMapPrimary, "Primary_hPercent"); | |
9f3c053c | 222 | |
223 | // Print full information on cells which are dead/noisy; | |
224 | // kTRUE -- print also the percentage % dead/% noisy | |
225 | PrintDeadNoisyCells(hBadCellMapPrimary, 0.9, 1.); // in 90-100% of runs | |
226 | // PrintDeadNoisyCells(hBadCellMapPrimary, 0.5, 0.9, kTRUE); // in 50-90% of runs | |
227 | // PrintDeadNoisyCells(hBadCellMapPrimary, 0.3, 0.5); // in 30-50% of runs | |
35ac8df7 | 228 | // PrintDeadNoisyCells(hBadCellMapPrimary, 0.1, 0.5); // in 10-50% of runs |
9f3c053c | 229 | |
230 | // visualize dead/noisy cell map for EMCAL/PHOS; requires aliroot | |
231 | DrawOccupancy(nruns, runNumbers, hBadCellMapPrimary, "hDeadNoisyCellsOccupancy"); | |
232 | ||
233 | // EMCAL: print full information on missing/noisy parts (e.g. RCUs); requires aliroot | |
234 | // PrintEMCALProblematicBlocks(nruns, runNumbers, hBadCellMapPrimary); | |
235 | ||
236 | ||
237 | /* RUNS QUALITY SECTION: CLUSTER AVERAGES + PI0 AVERAGES | |
238 | * | |
239 | */ | |
240 | ||
241 | // First, draw cluster averages per run; | |
242 | // 1 = minimum number of cells in a cluster; | |
243 | // 0.3GeV = minimum cluster energy; | |
244 | // -1 = maximum cluster energy = infinity (in fact, 20GeV in the default configuration). | |
245 | DrawClusterAveragesPerRun(nruns, runNumbers, 1, 0.3, -1); | |
246 | ||
247 | // Second, draw the same cluster averages per run, but corrected for detector acceptance | |
248 | DrawClusterAveragesPerRun(nruns, runNumbers, 1, 0.3, -1, hBadCellMapPrimary); | |
249 | ||
250 | // Draw random slices of the pi0 peak in one supermodule and in whole detector | |
35ac8df7 | 251 | // Int_t rndRun = runNumbers[gRandom->Integer(nruns)]; |
252 | Int_t rndRun = runNumbers[30]; | |
253 | DrawPi0Slice(rndRun, 1); | |
254 | DrawPi0Slice(rndRun, 2); | |
255 | DrawPi0Slice(rndRun, 3); | |
256 | DrawPi0Slice(rndRun, -1); | |
9f3c053c | 257 | |
258 | // Draw number of pi0s per event, pi0 mass and width | |
259 | DrawPi0Averages(nruns, runNumbers); | |
260 | ||
261 | // Draw pi0 values per event with pi0 photons both in the same supermodule | |
262 | // DrawPi0Averages(nruns, runNumbers, kTRUE); | |
263 | ||
264 | // Draw pi0 values per event with pi0 photons both in the same supermodule | |
265 | // + correct for supermodule acceptance (should not be treated to be 100% reliable!) | |
266 | // DrawPi0Averages(nruns, runNumbers, kTRUE, hBadCellMap[0]); | |
267 | ||
268 | // Draw pi0 values per event with pi0 | |
269 | // + correct for supermodule acceptance (should not be treated to be 100% reliable!) | |
c9365469 | 270 | // DrawPi0Averages(nruns, runNumbers, kFALSE, hBadCellMap[0]); |
9f3c053c | 271 | |
272 | // Draw pi0 distributions (helps to take decision on bad runs); | |
273 | // first argument -- suffix from canvas title; | |
274 | // second argument -- number of bins in distributions | |
275 | // DrawPi0Distributions("", 100); | |
276 | // DrawPi0Distributions("SM1", 100) | |
277 | // DrawPi0Distributions("_sameSM", 100); | |
278 | // DrawPi0Distributions("_sameSM_corr4accep", 100); | |
279 | ||
280 | // !!! | |
281 | return; | |
282 | ||
283 | /* SHAPE ANALYSIS | |
284 | * | |
285 | * Lazy/curious boundary is here! ------------------------- | |
286 | * Do not pass it if you do not fill curious enough! ;) | |
287 | * | |
288 | * The shape analysis below belongs to the main bad cell searching criteria. | |
289 | * However, it may require some extra work in order to make it usable. | |
290 | * | |
291 | * The idea behind is simple: fit the energy spectrum of a cell with | |
292 | * the function A*exp(-B*x)/x^2 (which proved to be a very good fit), | |
293 | * draw parameters A, B and chi2/ndf distributions and notice | |
294 | * cells which deviate too much from the averages. | |
295 | * | |
296 | * Huge statistics (>20M events) is necessary for this to work reliably. | |
297 | * | |
298 | * TODO: test on PHOS | |
299 | */ | |
300 | ||
301 | // The analysis below defines a cell as being bad simply if it is outside | |
302 | // of some good region for any of the parameters A, B, chi2/ndf. | |
303 | // The regions are depicted by vertical orange lines. The problem is that these | |
304 | // regions are usually not automatically determined correctly. | |
305 | // | |
306 | // The function syntax is the following: | |
307 | // 0.1-1.0 GeV -- fitting region; | |
308 | // hCellAmplitude -- the histogram which to take for processing, | |
309 | // the other possible choice is hCellAmplitudeNonLocMax; | |
310 | // 1000 -- number of bins is distributions. | |
311 | // | |
312 | // The next three groups of parameters are: | |
313 | // <text label> <maximum value on distribution> <left edge of the good region> <right edge of the good region> | |
314 | // | |
315 | // It is left/right edges which usually require manual tuning. | |
316 | // -1 means to initialize a parameter automatically. | |
317 | // | |
318 | TestCellShapes(0.25, 1., "hCellAmplitude", 1000, | |
319 | // maxval / left margin / right margin | |
320 | -1, -1,-1, // fit A | |
321 | -1, -1,-1, // fit B | |
322 | -1, -1,-1); // fit chi2/ndf | |
323 | ||
324 | ||
325 | /* DISTRIBUTION ANALYSIS | |
326 | * | |
327 | * Test for bad cells by plotting a distribution among cells. | |
328 | * | |
329 | * It is especially useful in searching for miscalibrated cells. | |
330 | * The function parameters are similar to parameters from shape analysis section. | |
331 | */ | |
332 | ||
333 | TestCellsTH1(nruns, runNumbers, "hCellLocMaxNTimesInClusterElow", | |
334 | "Number of times cell was local maximum, low energy, per event", "Entries", | |
335 | 400, -1,-1,-1); // nbins / maxval / left margin / right margin | |
336 | ||
337 | TestCellsTH1(nruns, runNumbers, "hCellLocMaxNTimesInClusterEhigh", | |
338 | "Number of times cell was local maximum, high energy, per event", "Entries", | |
339 | 400, -1,-1,-1); // nbins / maxval / left margin / right margin | |
340 | ||
341 | // TestCellsTH1(nruns, runNumbers, "hCellLocMaxETotalClusterElow", | |
342 | // "Total cluster energy for local maximum cell, low energy, per event", "Entries", | |
343 | // 400, -1,-1,-1); // nbins / maxval / left margin / right margin | |
344 | // | |
345 | // TestCellsTH1(nruns, runNumbers, "hCellLocMaxETotalClusterEhigh", | |
346 | // "Total cluster energy for local maximum cell, high energy, per event", "Entries", | |
347 | // 400, -1,-1,-1); // nbins / maxval / left margin / right margin | |
348 | ||
349 | // Three more tests for bad cells: | |
350 | // 1) total deposited energy; | |
351 | // 2) total number of entries; | |
352 | // 3) average energy = [total deposited energy]/[total number of entries]. | |
353 | // | |
354 | TestCellEandN(0.1, "hCellAmplitude", 1000, | |
355 | // maxval / left margin / right margin | |
356 | -1,-1,-1, // cell E | |
357 | -1,-1,-1, // cell N | |
358 | -1,-1,-1); // cell E/N | |
359 | ||
360 | // the same at high energies | |
361 | TestCellEandN(0.1, "hCellAmplitudeEhigh", 1000, | |
362 | // maxval / left margin / right margin | |
363 | -1,-1,-1, // cell E | |
364 | -1,-1,-1, // cell N | |
365 | -1,-1,-1); // cell E/N | |
366 | ||
367 | ||
368 | // The same as above, but for not a local maximum cells; require more statistics | |
369 | // TestCellsTH1(nruns, runNumbers, "hCellNonLocMaxNTimesInClusterElow", | |
370 | // "Number of times cell wasn't local maximum, low energy, per event", "Entries", | |
371 | // 400, -1,-1,-1); // nbins / maxval / left margin / right margin | |
372 | // | |
373 | // TestCellsTH1(nruns, runNumbers, "hCellNonLocMaxNTimesInClusterEhigh", | |
374 | // "Number of times cell wasn't local maximum, high energy, per event", "Entries", | |
375 | // 400, -1,-1,-1); // nbins / maxval / left margin / right margin | |
376 | // | |
377 | // TestCellsTH1(nruns, runNumbers, "hCellNonLocMaxETotalClusterElow", | |
378 | // "Total cluster energy for not local maximum cell, low energy, per event", "Entries", | |
379 | // 400, -1,-1,-1); // nbins / maxval / left margin / right margin | |
380 | // | |
381 | // TestCellsTH1(nruns, runNumbers, "hCellNonLocMaxETotalClusterEhigh", | |
382 | // "Total cluster energy for not local maximum cell, high energy, per event", "Entries", | |
383 | // 400, -1,-1,-1); // nbins / maxval / left margin / right margin | |
384 | // | |
385 | // TestCellEandN(0.1, "hCellAmplitudeNonLocMax", 1000, | |
386 | // // maxval / left margin / right margin | |
387 | // -1,-1,-1, // cell E | |
388 | // -1,-1,-1, // cell N | |
389 | // -1,-1,-1); // cell E/N | |
390 | // | |
391 | // TestCellEandN(0.1, "hCellAmplitudeEhighNonLocMax", 1000, | |
392 | // // maxval / left margin / right margin | |
393 | // -1,-1,-1, // cell E | |
394 | // -1,-1,-1, // cell N | |
395 | // -1,-1,-1); // cell E/N | |
396 | ||
397 | ||
398 | // TODO: cells time | |
399 | ||
400 | // The end ;) | |
401 | } | |
402 | ||
403 | ||
404 | ||
405 | /* BAD CELLS SEARCHING FUNCTIONS | |
406 | * | |
407 | */ | |
408 | ||
409 | //_________________________________________________________________________ | |
410 | TH2** FindDeadNoisyCellsPerRun(const Int_t nruns, Int_t runNumbers[], | |
411 | Double_t factorDead = 0.01, Double_t factorNoisy = 4., | |
412 | Int_t fnbins = 200, Double_t fxmax = 10.) | |
413 | { | |
414 | // Return bad cell candidate maps for the four criteria; | |
415 | // X axis -- cell absId, Y axis -- run index, content: 0=not bad, 1=noisy, -1=dead. | |
416 | // | |
417 | // For each run and each cell calculate factor = [cell value]/[averate over cells], | |
418 | // mark cell as noisy if factor >= factorNoisy, mark cell as dead if factor <= factorDead. | |
419 | // | |
420 | // fnbins, fxmax -- number of bins and X axis maximum value for factor distributions. | |
421 | // | |
422 | // Four criteria in order: | |
423 | char *hname[4] = {"hCellLocMaxNTimesInClusterElow", "hCellLocMaxNTimesInClusterEhigh", | |
424 | "hCellLocMaxETotalClusterElow", "hCellLocMaxETotalClusterEhigh"}; | |
425 | ||
426 | TH1* hFactorDistr[4]; | |
427 | TH2** hBadCellMap = new TH2*[4]; | |
428 | ||
429 | // take one histogram to get binning parameters | |
430 | TH1* one = (TH1*) gInputFile->Get(Form("run%i_%s",runNumbers[0],hname[0])); | |
431 | Int_t ncells = one->GetNbinsX(); | |
432 | Double_t amin = one->GetXaxis()->GetXmin(); | |
433 | Double_t amax = one->GetXaxis()->GetXmax(); | |
434 | ||
435 | // find dead/noisy cell candidates | |
436 | for (Int_t k = 0; k < 4; k++) { | |
437 | hBadCellMap[k] = new TH2C(Form("hBadCellMap_%s_fdead=%.3f_fnoisy=%.1f",hname[k],factorDead,factorNoisy), | |
438 | Form("Dead/noisy cell map (%s, fdead=%.3f, fnoisy=%.1f)",hname[k],factorDead,factorNoisy), | |
439 | ncells,amin,amax, nruns,0,nruns); | |
440 | hBadCellMap[k]->SetXTitle("AbsId"); | |
441 | hBadCellMap[k]->SetYTitle("Run index"); | |
442 | hBadCellMap[k]->SetTitleOffset(1.7,"Y"); | |
443 | ||
444 | hFactorDistr[k] = new TH1F(Form("hFactorDistr_%s_fdead=%.3f_fnoisy=%.1f", | |
445 | hname[k],factorDead,factorNoisy), "", fnbins,0,fxmax); | |
446 | hFactorDistr[k]->SetTitle("Factor distributions in all runs"); | |
447 | hFactorDistr[k]->SetXTitle("Factor"); | |
448 | hFactorDistr[k]->SetYTitle("Entries"); | |
449 | ||
450 | // run index loop | |
451 | for (Int_t ri = 0; ri < nruns; ri++) { | |
452 | TH1* one = (TH1*) gInputFile->Get(Form("run%i_%s", runNumbers[ri], hname[k])); | |
c9365469 | 453 | if (one == 0) { |
454 | Error("FindDeadNoisyCellsPerRun","Run %d does contain histogram %s",runNumbers[ri],hname[k]); | |
455 | continue; | |
456 | } | |
9f3c053c | 457 | // calculate average |
458 | Double_t av = 0; // average | |
459 | Int_t count = 0; // counted cells | |
460 | for (Int_t c = 1; c <= ncells; c++) { | |
461 | // do not count cells with zero content | |
462 | if (one->GetBinContent(c) == 0) continue; | |
463 | // do not count cells listed in the parameters section in the beginning of the macro | |
464 | if (IsCellMarked4Exclude(one->GetBinLowEdge(c))) continue; | |
465 | count++; | |
466 | av += one->GetBinContent(c); | |
467 | } | |
468 | ||
469 | // division by zero checks | |
470 | if (count == 0) { | |
471 | Warning("FindDeadNoisyCellsPerRun", Form("No cells counted at ri=%i",ri)); | |
472 | continue; | |
473 | } | |
474 | av /= count; | |
475 | ||
476 | if (av == 0) { | |
477 | Warning("FindDeadNoisyCellsPerRun", Form("Average is zero at ri=%i",ri)); | |
478 | continue; | |
479 | } | |
480 | ||
481 | // find dead/noisy candidates | |
482 | for (Int_t c = 1; c <= ncells; c++) { | |
483 | Double_t fac = one->GetBinContent(c)/av; | |
484 | hFactorDistr[k]->Fill(fac); | |
485 | ||
486 | if (fac <= factorDead) | |
487 | hBadCellMap[k]->SetBinContent(c, ri+1, -1); | |
488 | else if (fac >= factorNoisy) | |
489 | hBadCellMap[k]->SetBinContent(c, ri+1, 1); | |
490 | } | |
491 | ||
492 | delete one; | |
493 | } // run index loop | |
494 | } // criteria loop | |
495 | ||
496 | ||
497 | // draw factor distributions ... | |
498 | TCanvas *c1 = new TCanvas("hFactorDistr", "hFactorDistr", 400,400); | |
499 | gPad->SetLeftMargin(0.12); | |
500 | gPad->SetRightMargin(0.08); | |
501 | gPad->SetLogy(); | |
502 | hFactorDistr[0]->SetLineColor(kBlue); | |
503 | hFactorDistr[1]->SetLineColor(kRed); | |
504 | hFactorDistr[2]->SetLineColor(kGreen); | |
505 | hFactorDistr[3]->SetLineColor(kOrange); | |
506 | hFactorDistr[0]->Draw(); | |
507 | hFactorDistr[1]->Draw("same"); | |
508 | hFactorDistr[2]->Draw("same"); | |
509 | hFactorDistr[3]->Draw("same"); | |
510 | ||
511 | // ... with legend | |
512 | TLegend *leg = new TLegend(0.45,0.65,0.90,0.85); | |
513 | leg->AddEntry(hFactorDistr[0], "NTimesInCluster, low energy","l"); | |
514 | leg->AddEntry(hFactorDistr[2], "ETotalCluster, low energy","l"); | |
515 | leg->AddEntry(hFactorDistr[1], "NTimesInCluster, high energy","l"); | |
516 | leg->AddEntry(hFactorDistr[3], "EtotalCluster, high energy","l"); | |
517 | leg->Draw("same"); | |
518 | ||
519 | c1->Update(); | |
520 | ||
521 | return hBadCellMap; | |
522 | } | |
523 | ||
524 | //_________________________________________________________________________ | |
525 | void PrintCellsToExcludeFromAverages(TH2** hBadCellMap) | |
526 | { | |
527 | // Print cells suggested for exclusion from averages calculation | |
528 | ||
529 | Int_t ncells = hBadCellMap[0]->GetNbinsX(); | |
530 | Int_t nruns = hBadCellMap[0]->GetNbinsY(); | |
531 | ||
532 | Int_t *suggested = new Int_t[ncells]; | |
533 | memset(suggested, 0, ncells*sizeof(Int_t)); | |
534 | ||
535 | for (Int_t c = 1; c <= ncells; c++) | |
536 | for (Int_t ri = 1; ri <= nruns; ri++) | |
537 | if (hBadCellMap[0]->GetBinContent(c, ri) != 0 || hBadCellMap[2]->GetBinContent(c, ri) != 0 || | |
538 | hBadCellMap[1]->GetBinContent(c, ri) > 0 || hBadCellMap[3]->GetBinContent(c, ri) > 0) // NOTE: dead not counted | |
539 | suggested[c-1]++; | |
540 | ||
541 | printf("Suggested cells to switch off in averages calculations (copai 'n paste!):\n"); | |
542 | printf("Int_t excells[] = {"); | |
543 | Int_t n = 0; | |
544 | for (Int_t c = 1; c <= ncells; c++) | |
545 | if (suggested[c-1] >= 0.4*nruns) {// 40% of runs threshold | |
546 | printf("%s%i", n == 0 ? "" : ",", c); | |
547 | n++; | |
548 | } | |
549 | printf("};\nInt_t nexc = %i;\n\n",n); | |
550 | } | |
551 | ||
552 | //_________________________________________________________________________ | |
553 | Bool_t IsCellMarked4Exclude(Int_t absId) | |
554 | { | |
555 | // Return true if a cell is in excells[] array | |
556 | ||
557 | static TH1* hExclCells = NULL; | |
558 | ||
559 | // one time initialization | |
560 | if (!hExclCells) { | |
561 | hExclCells = new TH1F("hExclCells", "", 20000,0,20000); | |
562 | for (Int_t c = 0; c < nexc; c++) | |
563 | hExclCells->SetBinContent(hExclCells->FindBin(excells[c]), 1); | |
564 | } | |
565 | ||
566 | return (hExclCells->GetBinContent(hExclCells->FindBin(absId)) > 0 ? kTRUE : kFALSE); | |
567 | } | |
568 | ||
569 | //_________________________________________________________________________ | |
570 | void DrawDeadNoisyCellMap(TH2* hmap, char* cname) | |
571 | { | |
572 | // Visualize dead/noisy cell map; | |
573 | // cname -- canvas name. | |
574 | ||
c9365469 | 575 | // Define only 3 color: dead=blue, noisy=red, normal=white |
576 | Int_t color_index[]={kBlue,kWhite,kRed}; | |
577 | gStyle->SetPalette(3,color_index); | |
578 | ||
579 | TCanvas *c1 = new TCanvas(cname, cname, 0,0,1000,600); | |
580 | c1->Divide(1,3); | |
581 | TH2* hDummy[3]; | |
582 | ||
583 | // Use hardware numeration of modules (2,3,4) | |
584 | // instead of offline numeration (3,2,1) | |
585 | for (Int_t iM=1; iM<=3; iM++) { | |
586 | c1->cd(4-iM); | |
587 | gPad->SetLeftMargin(0.04); | |
588 | gPad->SetRightMargin(0.02); | |
589 | gPad->SetTopMargin(0.10); | |
590 | gPad->SetBottomMargin(0.13); | |
591 | hDummy[iM-1] = (TH2*) hmap->Clone(Form("hDummy_%s_mod%d",cname,iM)); | |
592 | hDummy[iM-1]->SetTitle(Form("Module %d",5-iM)); | |
593 | hDummy[iM-1]->SetLabelSize(0.08,"XY"); | |
594 | hDummy[iM-1]->SetTitleSize(0.08,"XY"); | |
595 | hDummy[iM-1]->SetTitleOffset(0.8,"X"); | |
596 | hDummy[iM-1]->SetTitleOffset(0.25,"Y"); | |
597 | hDummy[iM-1]->SetTickLength(0.01,"Y"); | |
598 | hDummy[iM-1]->SetNdivisions(520,"X"); | |
599 | hDummy[iM-1]->SetAxisRange(3584*(iM-1)+1, 3584*iM, "X"); | |
600 | hDummy[iM-1]->Draw("col"); | |
601 | } | |
9f3c053c | 602 | |
603 | c1->Update(); | |
604 | } | |
605 | ||
c9365469 | 606 | // //_________________________________________________________________________ |
607 | // void DrawDeadNoisyCellMap(TH2* hmap, char* cname) | |
608 | // { | |
609 | // // Visualize dead/noisy cell map; | |
610 | // // cname -- canvas name. | |
611 | ||
612 | // TCanvas *c1 = new TCanvas(cname, cname, 0,0,600,600); | |
613 | // gPad->SetLeftMargin(0.14); | |
614 | // gPad->SetRightMargin(0.06); | |
615 | ||
616 | // // draw dummy to initialize axis ranges | |
617 | // TH2* hDummy = (TH2*) hmap->Clone(Form("hDummy_%s",cname)); | |
618 | // hDummy->Reset(); | |
619 | // hDummy->Draw(); | |
620 | ||
621 | // for (Int_t c = 1; c <= hmap->GetNbinsX(); c++) //cell number | |
622 | // for (Int_t ri = 1; ri <= hmap->GetNbinsY(); ri++) { //run index | |
623 | // Double_t stat = hmap->GetBinContent(c, ri); // cell status | |
624 | // if (stat == 0) continue; | |
625 | ||
626 | // Double_t x = hmap->GetXaxis()->GetBinCenter(c); | |
627 | // Double_t y1 = hmap->GetYaxis()->GetBinLowEdge(ri); | |
628 | // Double_t y2 = hmap->GetYaxis()->GetBinUpEdge(ri); | |
629 | ||
630 | // // draw a line; FIXME: what is a better choice? | |
631 | // TLine* line = new TLine(x,y1,x,y2); | |
632 | // line->SetLineWidth(1); | |
633 | // if (stat > 0) line->SetLineColor(kRed); // noisy cell | |
634 | // else line->SetLineColor(kBlue); // dead cell | |
635 | // line->Draw(); | |
636 | // } | |
637 | ||
638 | // c1->Update(); | |
639 | // } | |
640 | ||
9f3c053c | 641 | //_________________________________________________________________________ |
642 | void DrawDeadNoisyCellPercent(Int_t nruns, Int_t runNumbers[], TH2* hmap, char* cname) | |
643 | { | |
644 | // Show percent of runs/events for each cell being problematic; | |
645 | // cname -- canvas name. | |
646 | ||
647 | // binning parameters | |
648 | Int_t ncells = hmap->GetNbinsX(); | |
649 | Double_t amin = hmap->GetXaxis()->GetXmin(); | |
650 | Double_t amax = hmap->GetXaxis()->GetXmax(); | |
651 | ||
652 | // number of times cell was dead/noisy; | |
653 | Int_t *ndeadRuns = new Int_t[ncells]; | |
654 | Int_t *nnoisyRuns = new Int_t[ncells]; | |
655 | Double_t *ndeadEvents = new Double_t[ncells]; | |
656 | Double_t *nnoisyEvents = new Double_t[ncells]; | |
657 | ||
658 | // fill arrays above | |
659 | for (Int_t c = 0; c < ncells; c++) { | |
660 | ndeadRuns[c] = 0; | |
661 | nnoisyRuns[c] = 0; | |
662 | ndeadEvents[c] = 0; | |
663 | nnoisyEvents[c] = 0; | |
664 | ||
665 | for (Int_t ri = 0; ri < nruns; ri++) { | |
666 | Double_t stat = hmap->GetBinContent(c+1, ri+1); // cell status | |
667 | Int_t nevents = GetNumberOfEvents(runNumbers[ri]); | |
668 | ||
669 | if (stat > 0) { | |
670 | nnoisyRuns[c]++; | |
671 | nnoisyEvents[c] += nevents; | |
672 | } | |
673 | else if (stat < 0) { | |
674 | ndeadRuns[c]++; | |
675 | ndeadEvents[c] += nevents; | |
676 | } | |
677 | } // run index loop | |
678 | } // cell loop | |
679 | ||
680 | // total number of events | |
681 | Double_t ntotal = GetTotalNumberOfEvents(nruns, runNumbers); | |
682 | ||
683 | TCanvas *c1 = new TCanvas(cname, cname, 0,0,600,600); | |
684 | gPad->SetLeftMargin(0.14); | |
685 | gPad->SetRightMargin(0.06); | |
686 | ||
687 | // draw dummy histogram to initialize canvas | |
688 | TH1* hDummy = new TH1F(Form("hDummy_%s",cname), hmap->GetTitle(), ncells,amin,amax); | |
689 | hDummy->SetAxisRange(0,100, "Y"); | |
690 | hDummy->SetXTitle("AbsId"); | |
691 | hDummy->SetYTitle("Percent"); | |
692 | hDummy->Draw(); | |
693 | ||
694 | // to fill legend | |
695 | TLine *line1 = NULL; | |
696 | TLine *line2 = NULL; | |
697 | TLine *line3 = NULL; | |
698 | TLine *line4 = NULL; | |
699 | ||
700 | // draw results, FIXME: is where a better way? | |
701 | for (Int_t c = 0; c < ncells; c++) { | |
702 | Double_t x = hmap->GetXaxis()->GetBinCenter(c+1); | |
703 | Double_t y1 = 100.*ndeadEvents[c]/ntotal; | |
704 | Double_t y2 = 100.*ndeadRuns[c]/nruns; | |
705 | ||
706 | // events, dead bar | |
707 | if (ndeadEvents[c] > 0) { | |
708 | line1 = new TLine(x, 0, x, y1); | |
709 | line1->SetLineWidth(1); | |
710 | line1->SetLineColor(kBlue); | |
711 | line1->Draw(); | |
712 | } | |
713 | ||
714 | // events, noisy bar | |
715 | if (nnoisyEvents[c] > 0) { | |
716 | line2 = new TLine(x, y1, x, y1 + 100.*nnoisyEvents[c]/ntotal); | |
717 | line2->SetLineWidth(1); | |
718 | line2->SetLineColor(kRed); | |
719 | line2->Draw(); | |
720 | } | |
721 | ||
722 | // runs, dead bar | |
723 | if (ndeadRuns[c] > 0) { | |
724 | line3 = new TLine(x, 0, x, y2); | |
725 | line3->SetLineWidth(1); | |
726 | line3->SetLineStyle(7); | |
727 | line3->SetLineColor(kBlack); | |
728 | line3->Draw(); | |
729 | } | |
730 | ||
731 | // runs, noisy bar | |
732 | if (nnoisyRuns[c] > 0) { | |
733 | line4 = new TLine(x, y2, x, y2 + 100.*nnoisyRuns[c]/nruns); | |
734 | line4->SetLineWidth(1); | |
735 | line4->SetLineStyle(7); | |
736 | line4->SetLineColor(kOrange); | |
737 | line4->Draw(); | |
738 | } | |
739 | } | |
740 | ||
741 | // legend | |
742 | TLegend *leg = new TLegend(0.65,0.65,0.92,0.75); | |
743 | if (line1) leg->AddEntry(line1, "% of events, dead","l"); | |
744 | if (line2) leg->AddEntry(line2, "% of events, noisy","l"); | |
745 | if (line3) leg->AddEntry(line3, "% of runs, dead","l"); | |
746 | if (line4) leg->AddEntry(line4, "% of runs, noisy","l"); | |
747 | leg->Draw("same"); | |
748 | ||
749 | c1->Update(); | |
750 | } | |
751 | ||
752 | //_________________________________________________________________________ | |
753 | TH1* ConvertMapRuns2Events(const Int_t nruns, Int_t runNumbers[], TH1* inhisto) | |
754 | { | |
755 | // Returns a histogram in which run index axis is converted to number of events | |
756 | // by making a variable axis bin width. | |
757 | // If inhisto inherits from TH2, convert Y axis; convert X axis otherwise. | |
758 | ||
759 | // bin widths | |
760 | Double_t *nevents = new Double_t[nruns+1]; | |
761 | nevents[0] = 0; | |
762 | for (Int_t ri = 0; ri < nruns; ri++) | |
763 | nevents[1+ri] = nevents[ri] + GetNumberOfEvents(runNumbers[ri]); | |
764 | ||
765 | TH1* outh = (TH1*) inhisto->Clone(Form("%s_Events",inhisto->GetName())); | |
766 | ||
767 | if (inhisto->InheritsFrom("TH2")) { | |
768 | outh->GetYaxis()->Set(nruns, nevents); | |
769 | outh->SetYTitle("Number of events"); | |
770 | outh->SetTitleOffset(1.7,"Y"); | |
771 | } | |
772 | else {// TH1 case | |
773 | outh->GetXaxis()->Set(nruns, nevents); | |
774 | outh->SetXTitle("Number of events"); | |
775 | } | |
776 | ||
777 | return outh; | |
778 | } | |
779 | ||
780 | //_________________________________________________________________________ | |
781 | TH2* CombineDeadNoisyElowEhigh(TH2* hmapElow, TH2* hmapEhigh) | |
782 | { | |
783 | // Combine two maps at low/high energy into one, | |
784 | // do not count dead map at high energy. | |
785 | // NOTE: if cell is dead at Elow and noisy and Ehigh, set content = 2. | |
786 | ||
787 | TH2* hmap_combined = (TH2*) hmapElow->Clone(Form("%s+%s",hmapElow->GetName(),hmapEhigh->GetName())); | |
788 | ||
789 | for (Int_t c = 1; c <= hmapElow->GetNbinsX(); c++) | |
790 | for (Int_t ri = 1; ri <= hmapElow->GetNbinsY(); ri++) | |
791 | if (hmapEhigh->GetBinContent(c, ri) > 0) { | |
792 | hmap_combined->SetBinContent(c, ri, 1); | |
793 | ||
794 | if (hmapElow->GetBinContent(c, ri) < 0) | |
795 | hmap_combined->SetBinContent(c, ri, 2); | |
796 | } | |
797 | ||
798 | hmap_combined->SetTitle("Dead/noisy cell map, combined"); | |
799 | ||
800 | return hmap_combined; | |
801 | } | |
802 | ||
803 | //_________________________________________________________________________ | |
804 | void PrintDeadNoisyCells(TH2* hmap, Double_t percent1 = 0.95, Double_t percent2 = 1., Bool_t kPrintPercentage = kFALSE) | |
805 | { | |
806 | // Print full information on dead/noisy cells which are present in | |
807 | // percent1-percent2 portion of runs (percent1 excluded, percent2 included). | |
808 | ||
809 | Int_t ncells = hmap->GetNbinsX(); | |
810 | Int_t nruns = hmap->GetNbinsY(); | |
811 | ||
812 | Int_t nbad = 0; | |
813 | ||
814 | printf("Dead/noisy cells in >%.1f%% and <=%.1f%% of runs:", 100*percent1, 100*percent2); | |
815 | ||
816 | for (Int_t c = 1; c <= ncells; c++) { | |
817 | Int_t nrdead = 0; // count number of runs for the current cell | |
818 | Int_t nrnoisy = 0; | |
819 | ||
820 | for (Int_t ri = 1; ri <= nruns; ri++) | |
821 | if (hmap->GetBinContent(c,ri) < 0) nrdead++; | |
822 | else if (hmap->GetBinContent(c,ri) > 0) nrnoisy++; | |
823 | ||
824 | if (nrdead+nrnoisy > percent1*nruns && nrdead+nrnoisy <= percent2*nruns) { | |
825 | printf(" %.0f", hmap->GetBinLowEdge(c)); | |
826 | if (kPrintPercentage) | |
827 | printf("(%-2.0f/%-2.0f)", 100.*nrdead/nruns, 100.*nrnoisy/nruns); | |
828 | nbad++; | |
829 | } | |
830 | } | |
831 | ||
832 | printf(" (%i total)\n\n", nbad); | |
833 | } | |
834 | ||
835 | //_________________________________________________________________________ | |
836 | void DrawOccupancy(Int_t nruns, Int_t runNumbers[], TH2* hmap, char* cname) | |
837 | { | |
838 | // Draw bad cell map for EMCAL or PHOS; | |
839 | // cname -- canvas name. | |
840 | ||
c9365469 | 841 | gStyle->SetPalette(1); |
842 | ||
9f3c053c | 843 | // guess detector |
844 | if (hmap->GetNbinsX() % 1152 == 0) | |
845 | DrawEMCALOccupancy(nruns, runNumbers, hmap, cname); | |
846 | else | |
847 | DrawPHOSOccupancy(nruns, runNumbers, hmap, cname); | |
848 | } | |
849 | ||
850 | //_________________________________________________________________________ | |
851 | void DrawEMCALOccupancy(Int_t nruns, Int_t runNumbers[], TH2* hmap, char* cname) | |
852 | { | |
853 | // Draw bad cell map for EMCAL; | |
854 | // cname -- canvas name. | |
855 | ||
856 | Int_t nmods = hmap->GetNbinsX()/1152; // number of supermodules | |
857 | Int_t vsize = ceil(nmods/2.); // vertical size in canvas | |
858 | Int_t lastSector = (nmods-1)/2; // for pad number calculation | |
859 | ||
860 | TCanvas *c1 = new TCanvas(cname, cname, 800,200*vsize); | |
861 | c1->Divide(2, vsize); | |
862 | ||
863 | for (Int_t sm = 0; sm < nmods; sm++) | |
864 | { | |
865 | // top left is SM0, bottom right is SM9 | |
866 | Int_t side = sm%2; | |
867 | Int_t sector = sm/2; | |
868 | c1->cd((lastSector - sector)*2 + side + 1); | |
869 | ||
870 | TH2* hSM = new TH2C(Form("hEMCALSM%i_%s",sm,cname), Form("SM%i",sm), 48,0,48, 24,0,24); | |
871 | hSM->SetXTitle("iEta"); | |
872 | hSM->SetYTitle("iPhi"); | |
873 | ||
874 | // loop over supermodule cells | |
875 | for (Int_t c = 0; c < 1152; c++) { | |
876 | Int_t absId = 1152 * sm + c; | |
877 | ||
878 | for (Int_t ri = 0; ri < nruns; ri++) { | |
879 | if (hmap->GetBinContent(hmap->FindBin(absId,ri)) == 0) continue; | |
880 | ||
881 | Int_t nModule, eta, phi; | |
882 | AbsId2EtaPhi(absId, nModule, eta, phi); | |
883 | hSM->Fill(eta,phi); | |
884 | } | |
885 | } | |
886 | ||
887 | hSM->Draw("colz"); | |
888 | } // supermodule loop | |
889 | ||
890 | c1->Update(); | |
891 | } | |
892 | ||
893 | //_________________________________________________________________________ | |
894 | void DrawPHOSOccupancy(Int_t nruns, Int_t runNumbers[], TH2* hmap, char* cname) | |
895 | { | |
896 | // Draw bad cell map for PHOS; | |
897 | // cname -- canvas name. | |
898 | ||
899 | Int_t nmods = hmap->GetNbinsX()/3584; // number of supermodules | |
900 | Int_t vsize = nmods; // vertical size in canvas | |
901 | ||
902 | TCanvas *c1 = new TCanvas(cname, cname, 64*10*vsize,56*10); | |
903 | c1->Divide(vsize,1); | |
904 | ||
905 | TFile *fBadMap = TFile::Open("PHOS_BadMap.root","recreate"); | |
906 | for (Int_t sm = 1; sm <= nmods; sm++) | |
907 | { | |
908 | c1->cd(sm); | |
909 | gPad->SetLeftMargin(0.10); | |
910 | gPad->SetRightMargin(0.15); | |
911 | gPad->SetTopMargin(0.05); | |
912 | gPad->SetBottomMargin(0.10); | |
913 | ||
c9365469 | 914 | TH2* hSM = new TH2C(Form("hPHOSSM%i_%s",sm,cname), Form("Module %i",5-sm), 64,0.5,64.5, 56,0.5,56.5); |
9f3c053c | 915 | hSM->SetXTitle("x_{cell}"); |
916 | hSM->SetYTitle("z_{cell} "); | |
917 | hSM->SetZTitle("N_{runs} "); | |
918 | ||
919 | // loop over supermodule cells | |
920 | for (Int_t c = 1; c <= 3584; c++) { | |
921 | Int_t absId = 3584*(sm-1) + c; | |
922 | ||
923 | for (Int_t ri = 0; ri < nruns; ri++) { | |
924 | if (hmap->GetBinContent(hmap->FindBin(absId,ri)) == 0) continue; | |
925 | ||
926 | Int_t nModule, xCell, zCell; | |
927 | AbsId2EtaPhi(absId, nModule, xCell, zCell, 1); | |
928 | hSM->Fill(xCell,zCell); | |
929 | } | |
930 | } | |
931 | hSM->Write(); | |
932 | hSM->DrawCopy("colz"); | |
933 | } // supermodule loop | |
934 | ||
935 | fBadMap->Close(); | |
936 | c1->Update(); | |
937 | } | |
938 | ||
939 | //_________________________________________________________________________ | |
940 | void PrintEMCALProblematicBlocks(Int_t nruns, Int_t runNumbers[], TH2* hmap) | |
941 | { | |
942 | // Print, on a per run basis, complete information about EMCAL missing | |
943 | // (or dead or noisy) blocks. Missing/noisy EMCAL atomic part is a 2x2 | |
944 | // block (288 blocks per supermodule). | |
945 | ||
946 | // number of supermodules | |
947 | Int_t nmods = hmap->GetNbinsX()/1152; | |
948 | ||
949 | Printf("Problematic (missing/dead/noisy) 2x2 block numbers in EMCAL:"); | |
950 | ||
951 | // run index loop | |
952 | for (Int_t ri = 0; ri < nruns; ri++) { | |
953 | Printf("Run %i:", runNumbers[ri]); | |
954 | ||
955 | // supermodule loop | |
956 | for (Int_t sm = 0; sm < nmods; sm++) { | |
957 | ||
958 | // will be filled with the number of missing cells (from 0 to 4) | |
959 | Int_t blk2x2[288]; | |
960 | memset(blk2x2, 0, 288*sizeof(Int_t)); | |
961 | ||
962 | for (Int_t c = 0; c < 1152; c++) { | |
963 | Int_t absId = 1152 * sm + c; | |
964 | ||
965 | // select problematic cells only | |
966 | if (hmap->GetBinContent(hmap->FindBin(absId,ri)) == 0) continue; | |
967 | ||
968 | Int_t nModule, eta, phi; | |
969 | AbsId2EtaPhi(c, nModule, eta, phi); | |
970 | blk2x2[nModule]++; | |
971 | } | |
972 | ||
973 | // calculate number of bad blocks | |
974 | Int_t nbad2x2 = 0; | |
975 | for (Int_t b = 0; b < 288; b++) | |
976 | if (blk2x2[b] == 4) nbad2x2++; | |
977 | ||
978 | // no bad | |
979 | if (nbad2x2 == 0) continue; | |
980 | ||
981 | printf(" SM%i:", sm); | |
982 | ||
983 | // whole supermodule | |
984 | if (nbad2x2 == 288) { | |
985 | printf(" missing the whole supermodule!\n"); | |
986 | continue; | |
987 | } | |
988 | ||
989 | // RCUs | |
990 | if (nbad2x2 >= 144) { | |
991 | Int_t nRCU[2]; | |
992 | nRCU[0] = 0; | |
993 | nRCU[1] = 0; | |
994 | ||
995 | for (Int_t b = 0; b < 288; b++) | |
996 | if (blk2x2[b] == 4) nRCU[GetEMCALRCUNumber(b)]++; | |
997 | ||
998 | if (nRCU[0] == 144) { | |
999 | printf(" RCU0"); | |
1000 | for (Int_t b = 0; b < 288; b++) | |
1001 | if (blk2x2[b] == 4 && GetEMCALRCUNumber(b) == 0) blk2x2[b] = 0; | |
1002 | } | |
1003 | ||
1004 | if (nRCU[1] == 144) { | |
1005 | printf(" RCU1"); | |
1006 | for (Int_t b = 0; b < 288; b++) | |
1007 | if (blk2x2[b] == 4 && GetEMCALRCUNumber(b) == 1) blk2x2[b] = 0; | |
1008 | } | |
1009 | ||
1010 | nbad2x2 -= 144; | |
1011 | } | |
1012 | ||
1013 | // the rest | |
1014 | if (nbad2x2 > 0) { | |
1015 | for (Int_t b = 0; b < 288; b++) | |
1016 | if (blk2x2[b] == 4) printf(" %i", b); | |
1017 | printf(" (%i)", nbad2x2); | |
1018 | } | |
1019 | ||
1020 | printf("\n"); | |
1021 | ||
1022 | } // supermodule loop | |
1023 | } // run index loop | |
1024 | ||
1025 | Printf(""); | |
1026 | } | |
1027 | ||
1028 | //_________________________________________________________________________ | |
1029 | Int_t GetEMCALRCUNumber(Int_t nModule) | |
1030 | { | |
1031 | // Returns RCU number for a 2x2 block in EMCAL; | |
1032 | // nModule -- block number (0-287). | |
1033 | ||
1034 | static TH1* hRCUs = NULL; | |
1035 | ||
1036 | // one-time initialization | |
1037 | if (!hRCUs) { | |
1038 | hRCUs = new TH1F("hRCU1", "", 288,0,288); | |
1039 | ||
1040 | Int_t RCU1[144] = {8,9,10,11,20,21,22,23,32,33,34,35,44,45,46,47,56,57,58,59,68,69,70,71,80,81,82,83, | |
1041 | 92,93,94,95,104,105,106,107,116,117,118,119,128,129,130,131,140,141,142,143,148,149,150,151, | |
1042 | 152,153,154,155,160,161,162,163,164,165,166,167,172,173,174,175,176,177,178,179,184,185,186, | |
1043 | 187,188,189,190,191,196,197,198,199,200,201,202,203,208,209,210,211,212,213,214,215,220,221, | |
1044 | 222,223,224,225,226,227,232,233,234,235,236,237,238,239,244,245,246,247,248,249,250,251,256, | |
1045 | 257,258,259,260,261,262,263,268,269,270,271,272,273,274,275,280,281,282,283,284,285,286,287}; | |
1046 | ||
1047 | for (Int_t i = 0; i < 144; i++) hRCUs->SetBinContent(RCU1[i]+1, 1); | |
1048 | } | |
1049 | ||
1050 | return hRCUs->GetBinContent(nModule+1); | |
1051 | } | |
1052 | ||
1053 | //_________________________________________________________________________ | |
1054 | void AbsId2EtaPhi(Int_t absId, Int_t &nModule, Int_t &eta, Int_t &phi, Int_t det = 0) | |
1055 | { | |
1056 | // Converts cell absId --> (sm,eta,phi); | |
1057 | // | |
1058 | // nModule -- 2x2 block number for EMCAL; | |
1059 | // det -- detector: 0/EMCAL, 1/PHOS. | |
1060 | ||
1061 | // EMCAL | |
1062 | if (det == 0) { | |
1063 | AliEMCALGeometry *geomEMCAL = AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1"); | |
1064 | ||
1065 | Int_t nSupMod, nIphi, nIeta; | |
1066 | geomEMCAL->GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta); | |
1067 | geomEMCAL->GetCellPhiEtaIndexInSModule(nSupMod, nModule, nIphi, nIeta, phi, eta); | |
1068 | return; | |
1069 | } | |
1070 | ||
1071 | // PHOS | |
1072 | else if (det == 1) { | |
1073 | AliPHOSGeometry *geomPHOS = AliPHOSGeometry::GetInstance("IHEP"); | |
1074 | ||
1075 | Int_t relid[4]; | |
1076 | geomPHOS->AbsToRelNumbering(absId, relid); | |
1077 | //sm = relid[0]; | |
1078 | eta = relid[2]; | |
1079 | phi = relid[3]; | |
1080 | return; | |
1081 | } | |
1082 | ||
1083 | // DCAL | |
1084 | // not implemented | |
1085 | // | |
1086 | ||
1087 | Error("AbsId2EtaPhi", "Wrong detector"); | |
1088 | abort(); | |
1089 | } | |
1090 | ||
1091 | //_________________________________________________________________________ | |
1092 | void TestCellsTH1(Int_t nruns, Int_t runNumbers[], char *hname, | |
1093 | char* htitle = "", char* hytitle = "", | |
1094 | Int_t dnbins = 200, Double_t dmaxval = -1, Double_t goodmin = -1, Double_t goodmax = -1) | |
1095 | { | |
1096 | // Test for bad cells by plotting a distribution of a TH1 histogram. | |
1097 | // The histogram is obtained as a sum over runs of TH1 per run. | |
1098 | // | |
1099 | // hname -- histogram name to process; | |
1100 | // htitle, hytitle -- histogram title and Y axis title; | |
1101 | // dnbins -- number of bins in distribution; | |
1102 | // dmaxval -- X axis maximum in distribution. | |
1103 | // goodmin,goodmax -- the region outside which a cell is considered as bad. | |
1104 | // | |
1105 | // -1 value for maxval/goodmin/goodmax -- process a variable automatically. | |
1106 | ||
1107 | // initialize histogram | |
1108 | TH1* histo = (TH1*) gInputFile->Get(Form("run%i_%s", runNumbers[0], hname)); | |
1109 | histo->SetName(hname); | |
1110 | histo->SetTitle(htitle); | |
1111 | histo->SetXTitle("AbsId"); | |
1112 | histo->SetYTitle(hytitle); | |
1113 | ||
1114 | // sum over runs | |
1115 | for (Int_t i = 1; i < nruns; i++) { | |
1116 | TH1* h = (TH1*) gInputFile->Get(Form("run%i_%s", runNumbers[i], hname)); | |
1117 | histo->Add(h); | |
1118 | delete h; | |
1119 | } | |
1120 | ||
1121 | histo->Scale(1./(Double_t)GetTotalNumberOfEvents(nruns, runNumbers)); | |
1122 | Process1(histo, Form("TestCellsTH1_%s",hname), dnbins, dmaxval, goodmin, goodmax); | |
1123 | } | |
1124 | ||
1125 | //_________________________________________________________________________ | |
1126 | void TestCellEandN(Double_t Emin = 0.1, char* hname = "hCellAmplitude", Int_t dnbins = 200, | |
1127 | Double_t maxval1 = -1, Double_t goodmin1 = -1, Double_t goodmax1 = -1, | |
1128 | Double_t maxval2 = -1, Double_t goodmin2 = -1, Double_t goodmax2 = -1, | |
1129 | Double_t maxval3 = -1, Double_t goodmin3 = -1, Double_t goodmax3 = -1) | |
1130 | { | |
1131 | // Three more tests for bad cells: | |
1132 | // 1) total deposited energy; | |
1133 | // 2) total number of entries; | |
1134 | // 3) average energy = [total deposited energy]/[total number of entries]. | |
1135 | // | |
1136 | // Based on summary histograms. Possible choises: | |
1137 | // hCellAmplitude, hCellAmplitudeEhigh, hCellAmplitudeNonLocMax, hCellAmplitudeEhighNonLocMax | |
1138 | // | |
1139 | // Emin -- minimum cell amplitude to count; | |
1140 | // hname -- name (in file) of TH2 histogram with cell amplitudes; | |
1141 | // dnbins -- number of bins in distributions; | |
1142 | // maxval[123] -- maximum values on distributions for the criteria 1),2),3) respectively; | |
1143 | // goodmin[123],goodmax[123] -- the regions on distributions outside those a cell is considered as bad. | |
1144 | ||
1145 | // input; X axis -- absId numbers | |
1146 | TH2 *hCellAmplitude = (TH2*) gInputFile->Get(hname); | |
1147 | ||
1148 | // binning parameters | |
1149 | Int_t ncells = hCellAmplitude->GetNbinsX(); | |
1150 | Double_t amin = hCellAmplitude->GetXaxis()->GetXmin(); | |
1151 | Double_t amax = hCellAmplitude->GetXaxis()->GetXmax(); | |
1152 | ||
1153 | TH1* hCellEtotal = new TH1F(Form("%s_hCellEtotal_E%.2f",hname,Emin), | |
1154 | Form("Total deposited energy, E > %.2f GeV",Emin), ncells,amin,amax); | |
1155 | hCellEtotal->SetXTitle("AbsId"); | |
1156 | hCellEtotal->SetYTitle("Energy, GeV"); | |
1157 | ||
1158 | TH1F *hCellNtotal = new TH1F(Form("%s_hCellNtotal_E%.2f",hname,Emin), | |
1159 | Form("Total number of entries, E > %.2f GeV",Emin), ncells,amin,amax); | |
1160 | hCellNtotal->SetXTitle("AbsId"); | |
1161 | hCellNtotal->SetYTitle("Entries"); | |
1162 | ||
1163 | TH1F *hCellEtoNtotal = new TH1F(Form("%s_hCellEtoNtotal_E%.2f",hname,Emin), | |
1164 | Form("Average energy per hit, E > %.2f GeV",Emin), ncells,amin,amax); | |
1165 | hCellEtoNtotal->SetXTitle("AbsId"); | |
1166 | hCellEtoNtotal->SetYTitle("Energy, GeV"); | |
1167 | ||
1168 | // fill cells | |
1169 | for (Int_t c = 1; c <= ncells; c++) { | |
1170 | Double_t Esum = 0; | |
1171 | Double_t Nsum = 0; | |
1172 | ||
1173 | for (Int_t i = 1; i <= hCellAmplitude->GetNbinsY(); i++) { | |
1174 | Double_t E = hCellAmplitude->GetYaxis()->GetBinCenter(i); | |
1175 | Double_t N = hCellAmplitude->GetBinContent(c, i); | |
1176 | if (E < Emin) continue; | |
1177 | Esum += E*N; | |
1178 | Nsum += N; | |
1179 | } | |
1180 | ||
1181 | hCellEtotal->SetBinContent(c, Esum); | |
1182 | hCellNtotal->SetBinContent(c, Nsum); | |
1183 | ||
1184 | if (Nsum > 0.5) // number of entries >= 1 | |
1185 | hCellEtoNtotal->SetBinContent(c, Esum/Nsum); | |
1186 | } | |
1187 | ||
1188 | delete hCellAmplitude; | |
1189 | ||
1190 | Process1(hCellEtotal, Form("%s_CellE", hname), dnbins, maxval1, goodmin1, goodmax1); | |
1191 | Process1(hCellNtotal, Form("%s_CellN", hname), dnbins, maxval2, goodmin2, goodmax2); | |
1192 | Process1(hCellEtoNtotal, Form("%s_CellE/N", hname), dnbins, maxval3, goodmin3, goodmax3); | |
1193 | } | |
1194 | ||
1195 | //_________________________________________________________________________ | |
1196 | void TestCellShapes(Double_t fitEmin, Double_t fitEmax, char* hname = "hCellAmplitude", Int_t dnbins = 1000, | |
1197 | Double_t maxval1 = -1, Double_t goodmin1 = -1, Double_t goodmax1 = -1, | |
1198 | Double_t maxval2 = -1, Double_t goodmin2 = -1, Double_t goodmax2 = -1, | |
1199 | Double_t maxval3 = -1, Double_t goodmin3 = -1, Double_t goodmax3 = -1) | |
1200 | { | |
1201 | // Test cells shape using fit function f(x)=A*exp(-B*x)/x^2. | |
1202 | // Produce values per cell + distributions for A, B and chi2/ndf parameters. | |
1203 | // | |
1204 | // fitEmin, fitEmax -- fit range; | |
1205 | // hname -- name (in file) of TH2 histogram with cell amplitudes; | |
1206 | // dnbins -- number of bins in distributions; | |
1207 | // maxval[123] -- maximum values on distributions for the criteria A, B and chi2/ndf respectively; | |
1208 | // goodmin[123],goodmax[123] -- the regions on distributions outside those a cell is considered as bad. | |
1209 | // | |
1210 | // -1 value for maxval/goodmin/goodmax -- process a variable automatically. | |
1211 | // | |
1212 | // Note: numbers are optimized for EMCAL. | |
1213 | // TODO: check for PHOS | |
1214 | ||
1215 | // input; X axis -- absId numbers | |
1216 | TH2 *hCellAmplitude = (TH2*) gInputFile->Get(hname); | |
1217 | ||
1218 | // binning parameters | |
1219 | Int_t ncells = hCellAmplitude->GetNbinsX(); | |
1220 | Double_t amin = hCellAmplitude->GetXaxis()->GetXmin(); | |
1221 | Double_t amax = hCellAmplitude->GetXaxis()->GetXmax(); | |
1222 | ||
1223 | // initialize histograms | |
1224 | TH1 *hFitA = new TH1F(Form("hFitA_%s",hname),"Fit A value", ncells,amin,amax); | |
1225 | hFitA->SetXTitle("AbsId"); | |
1226 | hFitA->SetYTitle("A"); | |
1227 | ||
1228 | TH1 *hFitB = new TH1F(Form("hFitB_%s",hname),"Fit B value", ncells,amin,amax); | |
1229 | hFitB->SetXTitle("AbsId"); | |
1230 | hFitB->SetYTitle("B"); | |
1231 | ||
1232 | TH1 *hFitChi2Ndf = new TH1F(Form("hFitChi2Ndf_%s",hname),"Fit #chi^{2}/ndf value", ncells,amin,amax); | |
1233 | hFitChi2Ndf->SetXTitle("AbsId"); | |
1234 | hFitChi2Ndf->SetYTitle("#chi^{2}/ndf"); | |
1235 | ||
1236 | // total number of events; to estimate A value | |
1237 | TH1* hNEventsProcessedPerRun = (TH1*) gInputFile->Get("hNEventsProcessedPerRun"); | |
1238 | Double_t ntotal = hNEventsProcessedPerRun->Integral(1, hNEventsProcessedPerRun->GetNbinsX()); | |
1239 | ||
1240 | // fitting function | |
1241 | TF1 *fit = new TF1("fit", "[0]*exp(-[1]*x)/x^2"); | |
1242 | fit->SetParLimits(0, ntotal*1e-8,ntotal*1e-4); | |
1243 | fit->SetParLimits(1, 0.,30.); | |
1244 | fit->SetParameter(0, ntotal*1e-6); | |
1245 | fit->SetParameter(1, 1.5); | |
1246 | ||
1247 | for (Int_t i = 1; i <= ncells; i++) { | |
1248 | TH1 *hCell = hCellAmplitude->ProjectionY("",i,i); | |
1249 | if (hCell->GetEntries() == 0) continue; | |
1250 | ||
1251 | hCell->Fit(fit, "0LQEM", "", fitEmin, fitEmax); | |
1252 | delete hCell; | |
1253 | ||
1254 | hFitA->SetBinContent(i, fit->GetParameter(0)); | |
1255 | hFitB->SetBinContent(i, fit->GetParameter(1)); | |
1256 | if (fit->GetNDF() != 0) | |
1257 | hFitChi2Ndf->SetBinContent(i, fit->GetChisquare()/fit->GetNDF()); | |
1258 | } | |
1259 | ||
1260 | delete hCellAmplitude; | |
1261 | ||
1262 | // automatic parameters, if requested | |
1263 | if (maxval1 < 0) maxval1 = 4e-6 * ntotal; | |
1264 | if (maxval2 < 0) maxval2 = 10.; | |
1265 | if (maxval3 < 0) maxval3 = 15.; | |
1266 | ||
1267 | Process1(hFitA, Form("%s_FitA", hname), dnbins, maxval1, goodmin1, goodmax1); | |
1268 | Process1(hFitB, Form("%s_FitB", hname), dnbins, maxval2, goodmin2, goodmax2); | |
1269 | Process1(hFitChi2Ndf, Form("%s_FitChi2ndf", hname), dnbins, maxval3, goodmin3, goodmax3); | |
1270 | } | |
1271 | ||
1272 | //_________________________________________________________________________ | |
1273 | void Process1(TH1* inhisto, char* label = "", Int_t dnbins = 200, | |
1274 | Double_t dmaxval = -1, Double_t goodmin = -1, Double_t goodmax = -1) | |
1275 | { | |
1276 | // Actual distribution analysis for a TH1 histogram: | |
1277 | // 1) create a distribution for the input histogram; | |
1278 | // 2) draw nicely; | |
1279 | // 3) take a decision about bad cells. | |
1280 | // | |
1281 | // inhisto -- input histogram; | |
1282 | // label -- text label for stdout; | |
1283 | // dnbins -- number of bins in distribution; | |
1284 | // goodmin,goodmax -- cells outside this region are considered as bad; | |
1285 | // dmaxval -- maximum value on distribution histogram. | |
1286 | // The later is required in cases where a bad cell kills the whole distribution: | |
1287 | // limiting distribution maximum value solves the problem. | |
1288 | ||
1289 | if (dmaxval < 0) | |
1290 | dmaxval = inhisto->GetMaximum()*1.01; // 1.01 - to see the last bin | |
1291 | ||
1292 | TH1 *distrib = new TH1F(Form("%sDistr",inhisto->GetName()), "", dnbins, inhisto->GetMinimum(), dmaxval); | |
1293 | distrib->SetXTitle(inhisto->GetYaxis()->GetTitle()); | |
1294 | distrib->SetYTitle("Entries"); | |
1295 | ||
1296 | // fill distribution | |
1297 | for (Int_t c = 1; c <= inhisto->GetNbinsX(); c++) | |
1298 | distrib->Fill(inhisto->GetBinContent(c)); | |
1299 | ||
1300 | // draw histogram + distribution | |
1301 | TCanvas *c1 = new TCanvas(inhisto->GetName(),inhisto->GetName(), 800,400); | |
1302 | c1->Divide(2,1); | |
1303 | ||
1304 | c1->cd(1); | |
1305 | gPad->SetLeftMargin(0.14); | |
1306 | gPad->SetRightMargin(0.06); | |
1307 | gPad->SetLogy(); | |
1308 | inhisto->SetTitleOffset(1.7,"Y"); | |
1309 | inhisto->Draw(); | |
1310 | ||
1311 | c1->cd(2); | |
1312 | gPad->SetLeftMargin(0.14); | |
1313 | gPad->SetRightMargin(0.10); | |
1314 | gPad->SetLogy(); | |
1315 | distrib->Draw(); | |
1316 | ||
1317 | // simple way to estimate the left margin for good cells region: | |
1318 | // go from left to right and find the first bin with content 2, | |
1319 | // then go from this bin right to left while bin content is nonzero | |
1320 | if (goodmin < 0) { | |
1321 | goodmin = distrib->GetXaxis()->GetXmin(); | |
1322 | ||
1323 | for (Int_t i = 1; i <= distrib->GetNbinsX(); i++) | |
1324 | if (distrib->GetBinContent(i) == 2) { | |
1325 | while (i > 1 && distrib->GetBinContent(i-1) > 0) i--; | |
1326 | goodmin = distrib->GetBinLowEdge(i); | |
1327 | break; | |
1328 | } | |
1329 | } | |
1330 | ||
1331 | // the same automatic algorithm as above, but reflected | |
1332 | if (goodmax < 0) { | |
1333 | goodmax = distrib->GetXaxis()->GetXmax(); | |
1334 | ||
1335 | for (Int_t i = distrib->GetNbinsX(); i >= 1; i--) | |
1336 | if (distrib->GetBinContent(i) == 2) { | |
1337 | while (i < distrib->GetNbinsX() && distrib->GetBinContent(i+1) > 0) i++; | |
1338 | goodmax = distrib->GetXaxis()->GetBinUpEdge(i); | |
1339 | break; | |
1340 | } | |
1341 | } | |
1342 | ||
1343 | // lines | |
1344 | TLine *lline = new TLine(goodmin, 0, goodmin, distrib->GetMaximum()); | |
1345 | lline->SetLineColor(kOrange); | |
1346 | lline->SetLineStyle(7); | |
1347 | lline->Draw(); | |
1348 | ||
1349 | TLine *rline = new TLine(goodmax, 0, goodmax, distrib->GetMaximum()); | |
1350 | rline->SetLineColor(kOrange); | |
1351 | rline->SetLineStyle(7); | |
1352 | rline->Draw(); | |
1353 | ||
1354 | // legend | |
1355 | TLegend *leg = new TLegend(0.60,0.82,0.9,0.88); | |
1356 | leg->AddEntry(lline, "Good region boundary","l"); | |
1357 | leg->Draw("same"); | |
1358 | ||
1359 | c1->Update(); | |
1360 | ||
1361 | ||
1362 | printf("Bad cells by criterum \"%s\":", label); | |
1363 | ||
1364 | // print bad cell numbers (outside goodmin,goodmax region) | |
1365 | Int_t ntot = 0; | |
1366 | for (Int_t c = 1; c <= inhisto->GetNbinsX(); c++) | |
1367 | if (inhisto->GetBinContent(c) < goodmin || inhisto->GetBinContent(c) > goodmax) { | |
1368 | printf(" %.0f", inhisto->GetBinLowEdge(c)); | |
1369 | ntot++; | |
1370 | } | |
1371 | ||
1372 | printf(" (%i total)\n\n", ntot); | |
1373 | } | |
1374 | ||
1375 | //_________________________________________________________________________ | |
1376 | void DrawCell(Int_t absId, Double_t fitEmin = 0.3, Double_t fitEmax = 1., | |
1377 | Double_t Emax = -1, char* hname = "hCellAmplitude") | |
1378 | { | |
1379 | // Draw one cell spectrum with a fit. | |
1380 | // | |
1381 | // fitEmin, fitEmax -- fit range; | |
1382 | // Emax -- maximum value on X axis to show (-1 = no limit); | |
1383 | // hname -- TH2 histogram name to process, possible values: | |
1384 | // "hCellAmplitude", "hCellAmplitudeEHigh", "hCellAmplitudeNonLocMax", "hCellAmplitudeEhighNonLocMax". | |
1385 | ||
1386 | // input; X axis -- absId numbers | |
1387 | TH2* hCellAmplitude = (TH2*) gInputFile->Get(hname); | |
1388 | ||
1389 | Int_t bin = hCellAmplitude->GetXaxis()->FindBin(absId); | |
1390 | TH1* hCell = hCellAmplitude->ProjectionY(Form("hCell%i_%s",absId,hname),bin,bin); | |
1391 | hCell->SetXTitle("Energy, GeV"); | |
1392 | hCell->SetYTitle("Entries"); | |
1393 | hCell->SetTitle(Form("Cell %i, %s", absId, hname)); | |
1394 | if (Emax > 0) hCell->SetAxisRange(0, Emax); | |
1395 | ||
1396 | // draw spectrum | |
1397 | TCanvas *c1 = new TCanvas(Form("hCell%i_%s",absId,hname), Form("hCell%i_%s",absId,hname), 400,400); | |
1398 | gPad->SetLeftMargin(0.12); | |
1399 | gPad->SetRightMargin(0.08); | |
1400 | gPad->SetBottomMargin(0.12); | |
1401 | gPad->SetLogy(); | |
1402 | hCell->Draw(); | |
1403 | ||
1404 | // total number of events | |
1405 | TH1* hNEventsProcessedPerRun = (TH1*) gInputFile->Get("hNEventsProcessedPerRun"); | |
1406 | Double_t ntotal = hNEventsProcessedPerRun->Integral(1, hNEventsProcessedPerRun->GetNbinsX()); | |
1407 | ||
1408 | // fit | |
1409 | TF1 *fit = new TF1("fit", "[0]*exp(-[1]*x)/x^2"); | |
1410 | fit->SetLineColor(kRed); | |
1411 | fit->SetLineWidth(2); | |
1412 | fit->SetParName(0, "A"); | |
1413 | fit->SetParName(1, "B"); | |
1414 | fit->SetParLimits(0, ntotal*1e-8,ntotal*1e-4); | |
1415 | fit->SetParLimits(1, 0.,30.); | |
1416 | fit->SetParameter(0, ntotal*1e-6); | |
1417 | fit->SetParameter(1, 1.); | |
1418 | hCell->Fit(fit,"LQEM", "", fitEmin, fitEmax); | |
1419 | ||
1420 | // legend | |
1421 | TLegend *leg = new TLegend(0.5,0.75,0.9,0.8); | |
1422 | leg->AddEntry(fit, "A*exp(-B*x)/x^{2}","l"); | |
1423 | leg->Draw("same"); | |
1424 | ||
1425 | c1->Update(); | |
1426 | } | |
1427 | ||
1428 | //_________________________________________________________________________ | |
1429 | void DrawCellTime(Int_t absId) | |
1430 | { | |
1431 | // Draw one cell time spectrum | |
1432 | ||
1433 | // input; X axis -- absId numbers | |
1434 | TH2* hCellTime = (TH2*) gInputFile->Get("hCellTime"); | |
1435 | ||
1436 | Int_t bin = hCellTime->GetXaxis()->FindBin(absId); | |
1437 | TH1* hCell = hCellTime->ProjectionY(Form("hCellTime%i",absId),bin,bin); | |
1438 | hCell->SetXTitle("Time, s"); | |
1439 | hCell->SetYTitle("Entries"); | |
1440 | hCell->SetTitle(Form("Cell %i time", absId)); | |
1441 | ||
1442 | // draw spectrum | |
1443 | TCanvas *c1 = new TCanvas(Form("hCellTime%i",absId), Form("hCellTime%i",absId), 400,400); | |
1444 | gPad->SetLeftMargin(0.12); | |
1445 | gPad->SetRightMargin(0.10); | |
1446 | gPad->SetBottomMargin(0.12); | |
1447 | gPad->SetLogy(); | |
1448 | hCell->Draw(); | |
1449 | ||
1450 | c1->Update(); | |
1451 | } | |
1452 | ||
1453 | ||
1454 | /* RUN AVERAGES AND RELATED FUNCTIONS | |
1455 | * | |
1456 | */ | |
1457 | ||
1458 | //_________________________________________________________________________ | |
1459 | void DrawClusterAveragesPerRun(Int_t nruns, Int_t runNumbers[], Int_t ncellsMin = 1, | |
1460 | Double_t eclusMin = 0.3, Double_t eclusMax = -1, | |
1461 | TH2* hmap = NULL) | |
1462 | { | |
1463 | // Draws cluster averages per run vs run index and number of events. | |
1464 | // NOTE: the averages are "smoothed" a little due to a finite bin width. | |
1465 | // | |
1466 | // ncellsMin -- minimum number of cells in cluster (>= 1); | |
1467 | // eclusMin,eclusMax -- minimum/maximum cluster energy cut | |
1468 | // (eclusMax = -1 means infinity); | |
1469 | // hmap -- dead/noisy cell map, which is used for acceptance calculation due | |
1470 | // to switched off detector parts. Acceptance is used to correct the average | |
1471 | // number of clusters per event. | |
1472 | ||
1473 | // names suffix | |
1474 | TString s(Form("_NC%i_Emin=%.2fGeV",ncellsMin,eclusMin)); | |
1475 | if (eclusMax > 0) s += TString(Form("_Emax=%.2fGeV",eclusMax)).Data(); | |
1476 | if (hmap) s += "_corr4accept"; | |
1477 | char *suff = s.Data(); | |
1478 | ||
1479 | if (eclusMax < 0) eclusMax = 1e+5; | |
1480 | ||
1481 | // supermodule region | |
1482 | Int_t SM1 = 0; | |
1483 | Int_t SM2 = 10; | |
1484 | while (SM1 <= SM2 && !gInputFile->Get(Form("run%i_hNCellsInClusterSM%i",runNumbers[0],SM1)) ) SM1++; | |
1485 | while (SM2 >= SM1 && !gInputFile->Get(Form("run%i_hNCellsInClusterSM%i",runNumbers[0],SM2)) ) SM2--; | |
1486 | ||
1487 | // initialize histograms | |
1488 | hAvECluster = new TH1F(Form("hAvECluster%s",suff), "Average cluster energy", nruns,0,nruns); | |
c9365469 | 1489 | SetRunLabel(hAvECluster,nruns,runNumbers,1); |
1490 | // hAvECluster->SetXTitle("Run index"); | |
1491 | hAvECluster->SetTitleOffset(0.3,"Y"); | |
9f3c053c | 1492 | hAvECluster->SetYTitle("Energy, GeV"); |
c9365469 | 1493 | hAvECluster->SetTickLength(0.01,"Y"); |
1494 | hAvECluster->SetLabelSize(0.06,"XY"); | |
1495 | hAvECluster->SetTitleSize(0.06,"XY"); | |
9f3c053c | 1496 | |
1497 | hAvNCluster = new TH1F(Form("hAvNCluster%s",suff), "Average number of clusters per event", nruns,0,nruns); | |
c9365469 | 1498 | SetRunLabel(hAvNCluster,nruns,runNumbers,1); |
1499 | // hAvNCluster->SetXTitle("Run index"); | |
1500 | hAvNCluster->SetTitleOffset(0.3,"Y"); | |
9f3c053c | 1501 | hAvNCluster->SetYTitle("Number of clusters"); |
c9365469 | 1502 | hAvNCluster->SetTickLength(0.01,"Y"); |
1503 | hAvNCluster->SetLabelSize(0.06,"XY"); | |
1504 | hAvNCluster->SetTitleSize(0.06,"XY"); | |
9f3c053c | 1505 | |
1506 | hAvNCellsInCluster = new TH1F(Form("hAvNCellsInCluster%s",suff), "Average number of cells in cluster", nruns,0,nruns); | |
c9365469 | 1507 | SetRunLabel(hAvNCellsInCluster,nruns,runNumbers,1); |
1508 | // hAvNCellsInCluster->SetXTitle("Run index"); | |
1509 | hAvNCellsInCluster->SetTitleOffset(0.3,"Y"); | |
9f3c053c | 1510 | hAvNCellsInCluster->SetYTitle("Number of cells"); |
c9365469 | 1511 | hAvNCellsInCluster->SetTickLength(0.01,"Y"); |
1512 | hAvNCellsInCluster->SetLabelSize(0.06,"XY"); | |
1513 | hAvNCellsInCluster->SetTitleSize(0.06,"XY"); | |
9f3c053c | 1514 | |
1515 | // initialize per SM histograms | |
1516 | TH1* hAvEClusterSM[10]; | |
1517 | TH1* hAvNClusterSM[10]; | |
1518 | TH1* hAvNCellsInClusterSM[10]; | |
1519 | ||
1520 | for (Int_t sm = SM1; sm <= SM2; sm++) { | |
1521 | hAvEClusterSM[sm] = new TH1F(Form("hAvEClusterSM%i%s",sm,suff),"", nruns,0,nruns); | |
1522 | hAvNClusterSM[sm] = new TH1F(Form("hAvNClusterSM%i%s",sm,suff),"", nruns,0,nruns); | |
1523 | hAvNCellsInClusterSM[sm] = new TH1F(Form("hAvNCellsInClusterSM%i%s",sm,suff),"", nruns,0,nruns); | |
1524 | } | |
1525 | ||
1526 | // fill all the histograms per run index | |
1527 | for (Int_t ri = 0; ri < nruns; ri++) | |
1528 | { | |
1529 | Int_t nevents = GetNumberOfEvents(runNumbers[ri]); | |
1530 | ||
1531 | // number of switched off supermodules | |
1532 | Int_t noSM = 0; | |
1533 | ||
1534 | Double_t Eclus_total = 0; // total cluster energy | |
1535 | Double_t Nclus_total = 0; // total number of clusters | |
1536 | Double_t Ncells_total = 0; // total number of cells | |
1537 | ||
1538 | // supermodule loop | |
1539 | for (Int_t sm = SM1; sm <= SM2; sm++) { | |
1540 | TH2* hNCellsInClusterSM = (TH2*) gInputFile->Get(Form("run%i_hNCellsInClusterSM%i",runNumbers[ri],sm)); | |
c9365469 | 1541 | if (hNCellsInClusterSM == 0) { |
1542 | Error("DrawClusterAveragesPerRun","Run %d does contain histogram %s",runNumbers[ri],Form("run%i_hNCellsInClusterSM%i",runNumbers[ri],sm)); | |
1543 | continue; | |
1544 | } | |
9f3c053c | 1545 | |
1546 | // the same as above, but per supermodule | |
1547 | Double_t Eclus_totalSM = 0; | |
1548 | Double_t Nclus_totalSM = 0; | |
1549 | Double_t Ncells_totalSM = 0; | |
1550 | ||
1551 | // X axis -- cluster energy, Y axis -- number of cells | |
1552 | for (Int_t x = 1; x <= hNCellsInClusterSM->GetNbinsX(); x++) | |
1553 | for (Int_t y = 1+ncellsMin; y <= hNCellsInClusterSM->GetNbinsY(); y++) {//NOTE: bin 1 correspond to ncellsMin=0 | |
1554 | Double_t Eclus = hNCellsInClusterSM->GetXaxis()->GetBinCenter(x); | |
1555 | Double_t Ncells = hNCellsInClusterSM->GetYaxis()->GetBinLowEdge(y); | |
1556 | Double_t Nclus = hNCellsInClusterSM->GetBinContent(x,y); | |
1557 | ||
1558 | // cut on cluster energy | |
1559 | if (Eclus < eclusMin || Eclus > eclusMax) continue; | |
1560 | ||
1561 | Eclus_totalSM += Eclus * Nclus; | |
1562 | Nclus_totalSM += Nclus; | |
1563 | Ncells_totalSM += Ncells * Nclus; | |
1564 | } | |
1565 | ||
1566 | delete hNCellsInClusterSM; | |
1567 | ||
1568 | // correct for acceptance | |
1569 | if (hmap) { | |
1570 | Double_t accep = GetAcceptance(sm, hmap, ri); | |
1571 | if (accep > 0) { | |
1572 | Eclus_totalSM /= accep; | |
1573 | Nclus_totalSM /= accep; | |
1574 | Ncells_totalSM /= accep; | |
1575 | } | |
1576 | else noSM++; | |
1577 | } | |
1578 | ||
1579 | Eclus_total += Eclus_totalSM; | |
1580 | Nclus_total += Nclus_totalSM; | |
1581 | Ncells_total += Ncells_totalSM; | |
1582 | ||
1583 | hAvNClusterSM[sm]->SetBinContent(ri+1, Nclus_totalSM/nevents); | |
1584 | if (Nclus_totalSM > 0) hAvEClusterSM[sm]->SetBinContent(ri+1, Eclus_totalSM/Nclus_totalSM); | |
1585 | if (Nclus_totalSM > 0) hAvNCellsInClusterSM[sm]->SetBinContent(ri+1, Ncells_totalSM/Nclus_totalSM); | |
1586 | } // supermodule loop | |
1587 | ||
1588 | hAvNCluster->SetBinContent(ri+1, Nclus_total/nevents/(SM2-SM1+1-noSM)); | |
1589 | if (Nclus_total > 0) hAvECluster->SetBinContent(ri+1, Eclus_total/Nclus_total); | |
1590 | if (Nclus_total > 0) hAvNCellsInCluster->SetBinContent(ri+1, Ncells_total/Nclus_total); | |
1591 | } // run index loop | |
1592 | ||
1593 | ||
1594 | /* Draw results vs run index | |
1595 | */ | |
1596 | ||
1597 | TCanvas *c1 = new TCanvas(Form("ClusterAveragesRuns%s",suff), | |
c9365469 | 1598 | Form("ClusterAveragesRuns%s",suff), 1000,707); |
1599 | c1->Divide(1,3); | |
9f3c053c | 1600 | |
1601 | // average cluster energy | |
1602 | c1->cd(1); | |
c9365469 | 1603 | gPad->SetLeftMargin(0.04); |
1604 | gPad->SetRightMargin(0.02); | |
1605 | gPad->SetTopMargin(0.10); | |
35ac8df7 | 1606 | gPad->SetBottomMargin(0.14); |
1607 | gPad->SetGridx(); | |
1608 | gPad->SetGridy(); | |
1609 | TLegend *leg = new TLegend(0.65,0.16,0.75,0.15+0.08*(SM2-SM1+1)); | |
9f3c053c | 1610 | |
1611 | hAvECluster->SetAxisRange(hAvECluster->GetMinimum()*0.5, hAvECluster->GetMaximum()*1.5,"Y"); | |
9f3c053c | 1612 | hAvECluster->SetLineWidth(2); |
1613 | hAvECluster->Draw(); | |
c9365469 | 1614 | leg->AddEntry(hAvECluster, "All Modules","l"); |
9f3c053c | 1615 | for (Int_t sm = SM1; sm <= SM2; sm++) { |
1616 | hAvEClusterSM[sm]->SetLineColor(colorsSM[sm]); | |
1617 | hAvEClusterSM[sm]->SetLineWidth(1); | |
1618 | hAvEClusterSM[sm]->Draw("same"); | |
c9365469 | 1619 | leg->AddEntry(hAvEClusterSM[sm], Form("Module %i",5-sm),"l"); |
9f3c053c | 1620 | } |
1621 | hAvECluster->Draw("same"); // to top | |
1622 | leg->Draw("same"); | |
1623 | ||
1624 | // average number of clusters | |
1625 | c1->cd(2); | |
c9365469 | 1626 | gPad->SetLeftMargin(0.04); |
1627 | gPad->SetRightMargin(0.02); | |
1628 | gPad->SetTopMargin(0.10); | |
35ac8df7 | 1629 | gPad->SetBottomMargin(0.14); |
1630 | gPad->SetGridx(); | |
1631 | gPad->SetGridy(); | |
1632 | leg = new TLegend(0.65,0.16,0.75,0.15+0.08*(SM2-SM1+1)); | |
9f3c053c | 1633 | |
1634 | hAvNCluster->SetAxisRange(0, hAvNCluster->GetMaximum()*1.5,"Y"); | |
9f3c053c | 1635 | hAvNCluster->SetLineWidth(2); |
1636 | hAvNCluster->Draw(); | |
c9365469 | 1637 | leg->AddEntry(hAvNCluster, Form("(All Moduls)/%i",SM2-SM1+1),"l"); |
9f3c053c | 1638 | for (Int_t sm = SM1; sm <= SM2; sm++) { |
1639 | hAvNClusterSM[sm]->SetLineColor(colorsSM[sm]); | |
1640 | hAvNClusterSM[sm]->SetLineWidth(1); | |
1641 | hAvNClusterSM[sm]->Draw("same"); | |
c9365469 | 1642 | leg->AddEntry(hAvNClusterSM[sm], Form("Module %i",5-sm),"l"); |
9f3c053c | 1643 | } |
1644 | hAvNCluster->Draw("same"); // to top | |
1645 | leg->Draw("same"); | |
1646 | ||
1647 | // average number of cells in cluster | |
1648 | c1->cd(3); | |
c9365469 | 1649 | gPad->SetLeftMargin(0.04); |
1650 | gPad->SetRightMargin(0.02); | |
1651 | gPad->SetTopMargin(0.10); | |
35ac8df7 | 1652 | gPad->SetBottomMargin(0.14); |
1653 | gPad->SetGridx(); | |
1654 | gPad->SetGridy(); | |
1655 | leg = new TLegend(0.65,0.16,0.75,0.15+0.08*(SM2-SM1+1)); | |
9f3c053c | 1656 | |
1657 | hAvNCellsInCluster->SetAxisRange(0, hAvNCellsInCluster->GetMaximum()*1.3,"Y"); | |
9f3c053c | 1658 | hAvNCellsInCluster->SetLineWidth(2); |
1659 | hAvNCellsInCluster->Draw(); | |
c9365469 | 1660 | leg->AddEntry(hAvNCellsInCluster, "All Modules","l"); |
9f3c053c | 1661 | for (Int_t sm = SM1; sm <= SM2; sm++) { |
1662 | hAvNCellsInClusterSM[sm]->SetLineColor(colorsSM[sm]); | |
1663 | hAvNCellsInClusterSM[sm]->SetLineWidth(1); | |
1664 | hAvNCellsInClusterSM[sm]->Draw("same"); | |
c9365469 | 1665 | leg->AddEntry(hAvNCellsInClusterSM[sm], Form("Module %i",5-sm),"l"); |
9f3c053c | 1666 | } |
1667 | hAvNCellsInCluster->Draw("same"); // to top | |
1668 | leg->Draw("same"); | |
35ac8df7 | 1669 | c1->Update(); |
9f3c053c | 1670 | |
35ac8df7 | 1671 | if (!trendPlotEvents) return; |
9f3c053c | 1672 | |
1673 | /* Draw the same vs number of events | |
1674 | */ | |
1675 | ||
1676 | TCanvas *c2 = new TCanvas(Form("ClusterAveragesEvents%s",suff), | |
c9365469 | 1677 | Form("ClusterAveragesEvents%s",suff), 1000,707); |
1678 | c2->Divide(1,3); | |
9f3c053c | 1679 | |
1680 | // average cluster energy | |
1681 | c2->cd(1); | |
c9365469 | 1682 | gPad->SetLeftMargin(0.04); |
1683 | gPad->SetRightMargin(0.02); | |
1684 | gPad->SetTopMargin(0.10); | |
1685 | gPad->SetBottomMargin(0.13); | |
9f3c053c | 1686 | leg = new TLegend(0.65,0.15,0.85,0.15+0.04*(SM2-SM1+1)); |
1687 | ||
1688 | TH1* hAvEClusterEvents = ConvertMapRuns2Events(nruns, runNumbers, hAvECluster); | |
1689 | hAvEClusterEvents->Draw(); | |
c9365469 | 1690 | leg->AddEntry(hAvEClusterEvents, "All Modules","l"); |
9f3c053c | 1691 | for (Int_t sm = SM1; sm <= SM2; sm++) { |
1692 | TH1* hAvEClusterSMEvents = ConvertMapRuns2Events(nruns, runNumbers, hAvEClusterSM[sm]); | |
1693 | hAvEClusterSMEvents->Draw("same"); | |
c9365469 | 1694 | leg->AddEntry(hAvEClusterSMEvents, Form("Module %i",5-sm),"l"); |
9f3c053c | 1695 | } |
1696 | hAvEClusterEvents->Draw("same"); // to top | |
1697 | leg->Draw("same"); | |
1698 | ||
1699 | // average number of clusters | |
1700 | c2->cd(2); | |
c9365469 | 1701 | gPad->SetLeftMargin(0.04); |
1702 | gPad->SetRightMargin(0.02); | |
1703 | gPad->SetTopMargin(0.10); | |
1704 | gPad->SetBottomMargin(0.13); | |
9f3c053c | 1705 | leg = new TLegend(0.65,0.15,0.85,0.15+0.04*(SM2-SM1+1)); |
1706 | ||
1707 | TH1* hAvNClusterEvents = ConvertMapRuns2Events(nruns, runNumbers, hAvNCluster); | |
1708 | hAvNClusterEvents->Draw(); | |
c9365469 | 1709 | leg->AddEntry(hAvNClusterEvents, Form("(All Modules)/%i",SM2-SM1+1),"l"); |
9f3c053c | 1710 | for (Int_t sm = SM1; sm <= SM2; sm++) { |
1711 | TH1* hAvNClusterSMEvents = ConvertMapRuns2Events(nruns, runNumbers, hAvNClusterSM[sm]); | |
1712 | hAvNClusterSMEvents->Draw("same"); | |
c9365469 | 1713 | leg->AddEntry(hAvNClusterSMEvents, Form("Module %i",5-sm),"l"); |
9f3c053c | 1714 | } |
1715 | hAvNClusterEvents->Draw("same"); // to top | |
1716 | leg->Draw("same"); | |
1717 | ||
1718 | // average number of cells in cluster | |
1719 | c2->cd(3); | |
c9365469 | 1720 | gPad->SetLeftMargin(0.04); |
1721 | gPad->SetRightMargin(0.02); | |
1722 | gPad->SetTopMargin(0.10); | |
1723 | gPad->SetBottomMargin(0.13); | |
9f3c053c | 1724 | leg = new TLegend(0.65,0.15,0.85,0.15+0.04*(SM2-SM1+1)); |
1725 | ||
1726 | TH1* hAvNCellsInClusterEvents = ConvertMapRuns2Events(nruns, runNumbers, hAvNCellsInCluster); | |
1727 | hAvNCellsInClusterEvents->Draw(); | |
c9365469 | 1728 | leg->AddEntry(hAvNCellsInClusterEvents, "All Modules","l"); |
9f3c053c | 1729 | for (Int_t sm = SM1; sm <= SM2; sm++) { |
1730 | TH1* hAvNCellsInClusterSMEvents = ConvertMapRuns2Events(nruns, runNumbers, hAvNCellsInClusterSM[sm]); | |
1731 | hAvNCellsInClusterSMEvents->Draw("same"); | |
c9365469 | 1732 | leg->AddEntry(hAvNCellsInClusterSMEvents, Form("Module %i",5-sm),"l"); |
9f3c053c | 1733 | } |
1734 | hAvNCellsInClusterEvents->Draw("same"); // to top | |
1735 | leg->Draw("same"); | |
1736 | ||
9f3c053c | 1737 | c2->Update(); |
1738 | } | |
1739 | ||
1740 | //_________________________________________________________________________ | |
1741 | Double_t GetAcceptance(Int_t sm, TH2* hmap, Int_t ri) | |
1742 | { | |
1743 | // Returns [#cells - #dead]/#cells for a supermodule. | |
1744 | // hmap -- dead/noisy cell map; | |
1745 | // ri -- run index. | |
1746 | ||
1747 | // guess number of cells per supermodule | |
1748 | Int_t nSM = 1152; // EMCAL | |
1749 | if (hmap->GetXaxis()->GetXmin() == 1) {// PHOS | |
1750 | nSM = 3584; | |
1751 | sm--; // starts from 1, convenient from 0 | |
1752 | } | |
1753 | ||
1754 | // count dead cells | |
1755 | Int_t ndead = 0; | |
1756 | for (Int_t k = 1; k <= nSM; k++) | |
1757 | if (hmap->GetBinContent(nSM*sm + k, ri+1) < 0) | |
1758 | ndead++; | |
1759 | ||
1760 | return ((Double_t) (nSM - ndead))/nSM; | |
1761 | } | |
1762 | ||
1763 | //_________________________________________________________________________ | |
1764 | void DrawPi0Slice(Int_t run, Int_t sm = -1) | |
1765 | { | |
1766 | // Draw the pi0 peak; | |
1767 | // run,sm -- run number and supermodule to take; | |
1768 | // sm < 0 -- draw for whole detector. | |
1769 | ||
1770 | TH1* h = NULL; | |
1771 | if (sm >= 0) {//particular supermodule | |
1772 | h = (TH1*) gInputFile->Get(Form("run%i_hPi0MassSM%iSM%i",run,sm,sm)); | |
c9365469 | 1773 | if (h == 0) { |
1774 | Error("DrawPi0Slice","Run %d does contain histogram %s",run,Form("run%i_hPi0MassSM%iSM%i",run,sm,sm)); | |
1775 | continue; | |
1776 | } | |
9f3c053c | 1777 | h->SetName(Form("hPi0SliceSM%i_run%i",sm,run)); |
35ac8df7 | 1778 | h->SetTitle(Form("#pi^{0} in M%i, run %i, %.0fk events", 5-sm, run, GetNumberOfEvents(run)/1e+3)); |
9f3c053c | 1779 | } |
1780 | else {// whole detector | |
1781 | for (Int_t sm1 = 0; sm1 < 10; sm1++) | |
1782 | for (Int_t sm2 = sm1; sm2 < 10; sm2++) { | |
1783 | TH1* one = (TH1*) gInputFile->Get(Form("run%i_hPi0MassSM%iSM%i",run,sm1,sm2)); | |
1784 | if (!one) continue; | |
1785 | ||
1786 | if (!h) h = one; | |
1787 | else { | |
1788 | h->Add(one); | |
1789 | delete one; | |
1790 | } | |
1791 | } | |
1792 | h->SetName(Form("hPi0Slice_run%i",run)); | |
1793 | h->SetTitle(Form("#pi^{0} in all modules, run %i, %.0fk events", run, GetNumberOfEvents(run)/1e+3)); | |
1794 | } | |
1795 | ||
1796 | h->SetXTitle("M_{#gamma#gamma}, GeV"); | |
1797 | h->SetYTitle("Entries"); | |
1798 | h->SetTitleOffset(1.7,"Y"); | |
1799 | ||
1800 | TCanvas *c1; | |
1801 | if (sm >= 0) c1 = new TCanvas(Form("hPi0SliceSM%i_run%i",sm,run),Form("hPi0SliceSM%i_run%i",sm,run), 400,400); | |
1802 | else c1 = new TCanvas(Form("hPi0Slice_run%i",run),Form("hPi0Slice_run%i",run), 400,400); | |
1803 | ||
1804 | gPad->SetLeftMargin(0.14); | |
1805 | gPad->SetRightMargin(0.06); | |
1806 | h->Draw(); | |
1807 | ||
1808 | Double_t nraw, enraw, mass, emass, sigma, esigma; | |
1809 | FitPi0(h, nraw, enraw, mass, emass, sigma, esigma); | |
1810 | ||
1811 | // draw background | |
1812 | TF1* fitfun = h->GetFunction("fitfun"); | |
1813 | if (fitfun) { | |
1814 | Double_t emin, emax; | |
1815 | fitfun->GetRange(emin, emax); | |
1816 | ||
1817 | backgr = new TF1("mypol2", "[0] + [1]*(x-0.135) + [2]*(x-0.135)^2", emin, emax); | |
1818 | backgr->SetLineColor(kBlue); | |
1819 | backgr->SetLineWidth(2); | |
1820 | backgr->SetLineStyle(3); | |
1821 | backgr->SetParameters(fitfun->GetParameter(3), fitfun->GetParameter(4), fitfun->GetParameter(5)); | |
1822 | backgr->Draw("same"); | |
1823 | } | |
1824 | ||
1825 | c1->Update(); | |
1826 | } | |
1827 | ||
1828 | //_________________________________________________________________________ | |
1829 | void DrawPi0Averages(Int_t nruns, Int_t runNumbers[], Bool_t samesm = kFALSE, TH2* hmap = NULL) | |
1830 | { | |
1831 | // Draw average number of pi0s per event, pi0 mass position and pi0 width per run index. | |
1832 | // Errors are also drawn. | |
1833 | // | |
1834 | // samesm -- take only pi0s for which gammas were is same SM; | |
1835 | // hmap -- if not NULL, do simple (area law based) acceptance correction. | |
1836 | // | |
1837 | // TODO: PHOS needs pi0 between SMs rather than in one SM | |
1838 | ||
1839 | // suffix to names | |
1840 | char *suff = TString(Form("%s%s", samesm ? "_sameSM":"", hmap ? "_corr4accep":"")).Data(); | |
1841 | ||
1842 | // supermodule region | |
1843 | Int_t SM1 = 0; | |
1844 | Int_t SM2 = 10; | |
1845 | while (SM1 <= SM2 && !gInputFile->Get(Form("run%i_hPi0MassSM%iSM%i",runNumbers[0],SM1,SM1)) ) SM1++; | |
1846 | while (SM2 >= SM1 && !gInputFile->Get(Form("run%i_hPi0MassSM%iSM%i",runNumbers[0],SM2,SM2)) ) SM2--; | |
1847 | ||
1848 | // initialize histograms for the entire detector | |
1849 | TH1* hPi0Num = new TH1F(Form("hPi0Num%s",suff),"Average number of #pi^{0}s per event", nruns,0,nruns); | |
c9365469 | 1850 | SetRunLabel(hPi0Num,nruns,runNumbers,1); |
1851 | // hPi0Num->SetXTitle("Run index"); | |
9f3c053c | 1852 | hPi0Num->SetYTitle("Number of #pi^{0}s"); |
1853 | ||
1854 | TH1* hPi0Mass = new TH1F(Form("hPi0Mass%s",suff),"#pi^{0} mass position", nruns,0,nruns); | |
c9365469 | 1855 | SetRunLabel(hPi0Mass,nruns,runNumbers,1); |
1856 | // hPi0Mass->SetXTitle("Run index"); | |
35ac8df7 | 1857 | hPi0Mass->SetYTitle("M_{#pi^{0}}, GeV/c^{2}"); |
9f3c053c | 1858 | |
1859 | TH1* hPi0Sigma = new TH1F(Form("hPi0Sigma%s",suff),"#pi^{0} width", nruns,0,nruns); | |
c9365469 | 1860 | SetRunLabel(hPi0Sigma,nruns,runNumbers,1); |
1861 | // hPi0Sigma->SetXTitle("Run index"); | |
35ac8df7 | 1862 | hPi0Sigma->SetYTitle("#sigma_{#pi^{0}}, GeV/c^{2}"); |
9f3c053c | 1863 | |
1864 | // initialize histograms per SM | |
1865 | TH1* hPi0NumSM[10]; | |
1866 | TH1* hPi0MassSM[10]; | |
1867 | TH1* hPi0SigmaSM[10]; | |
1868 | ||
1869 | for (Int_t sm = SM1; sm <= SM2; sm++) { | |
35ac8df7 | 1870 | hPi0NumSM[sm] = new TH1F(Form("hPi0NumSM%i%s",sm,suff),"", nruns,0,nruns); |
1871 | hPi0MassSM[sm] = new TH1F(Form("hPi0MassSM%i%s",sm,suff),"", nruns,0,nruns); | |
9f3c053c | 1872 | hPi0SigmaSM[sm] = new TH1F(Form("hPi0SigmaSM%i%s",sm,suff),"", nruns,0,nruns); |
1873 | } | |
1874 | ||
35ac8df7 | 1875 | TCanvas cFit("cFit","cFit"); |
9f3c053c | 1876 | // run index loop |
1877 | for (Int_t ri = 0; ri < nruns; ri++) | |
1878 | { | |
35ac8df7 | 1879 | fprintf(stderr, "\rDrawPi0Averages(): analysing run index %i/%i, run # %i... ", ri, nruns-1,runNumbers[ri]); |
9f3c053c | 1880 | |
1881 | Int_t nevents = GetNumberOfEvents(runNumbers[ri]); | |
1882 | ||
1883 | // per SM histos | |
1884 | for (Int_t sm = SM1; sm <= SM2; sm++) { | |
1885 | TH1* h = (TH1*) gInputFile->Get(Form("run%i_hPi0MassSM%iSM%i",runNumbers[ri],sm,sm)); | |
c9365469 | 1886 | if (h == 0) { |
1887 | Error("DrawPi0Averages","Run %d does contain histogram %s",runNumbers[ri],Form("run%i_hPi0MassSM%iSM%i",runNumbers[ri],sm,sm)); | |
1888 | continue; | |
1889 | } | |
9f3c053c | 1890 | |
1891 | // supermodule acceptance | |
1892 | Double_t accep = 1.; | |
1893 | if (hmap) accep = GetAcceptance(sm, hmap, ri); | |
1894 | if (accep == 0) continue; // missing SM | |
1895 | ||
1896 | Double_t nraw, enraw, mass, emass, sigma, esigma; | |
1897 | FitPi0(h, nraw, enraw, mass, emass, sigma, esigma); | |
1898 | ||
35ac8df7 | 1899 | hPi0NumSM[sm] ->SetBinContent(ri+1, nraw/nevents/accep); |
1900 | hPi0MassSM[sm] ->SetBinContent(ri+1, mass); | |
9f3c053c | 1901 | hPi0SigmaSM[sm]->SetBinContent(ri+1, sigma); |
1902 | ||
35ac8df7 | 1903 | hPi0NumSM[sm] ->SetBinError(ri+1, enraw/nevents/accep); |
1904 | hPi0MassSM[sm] ->SetBinError(ri+1, emass); | |
9f3c053c | 1905 | hPi0SigmaSM[sm]->SetBinError(ri+1, esigma); |
1906 | ||
1907 | delete h; | |
1908 | } // supermodule loop | |
1909 | ||
1910 | ||
1911 | /* fill for the entire detector | |
1912 | */ | |
1913 | TH1* hsum = (TH1*) gInputFile->Get(Form("run%i_hPi0MassSM%iSM%i",runNumbers[ri],SM1,SM1)); | |
c9365469 | 1914 | if (hsum == 0) { |
1915 | Error("DrawPi0Averages","Run %d does contain histogram %s",runNumbers[ri],Form("run%i_hPi0MassSM%iSM%i",runNumbers[ri],SM1,SM1)); | |
1916 | continue; | |
1917 | } | |
9f3c053c | 1918 | hsum->SetName("hSumTMP"); |
1919 | ||
1920 | // for acceptance correction | |
1921 | Int_t noSM = 0; | |
1922 | ||
1923 | for (Int_t sm1 = SM1; sm1 <= SM2; sm1++) | |
1924 | for (Int_t sm2 = sm1; sm2 <= SM2; sm2++) { | |
1925 | if (sm1 == SM1 && sm2 == SM2) continue; | |
1926 | if (samesm && sm1 != sm2) continue; | |
1927 | ||
1928 | TH1* h = (TH1*) gInputFile->Get(Form("run%i_hPi0MassSM%iSM%i",runNumbers[ri],sm1,sm2)); | |
c9365469 | 1929 | if (h == 0) { |
1930 | Error("DrawPi0Averages","Run %d does contain histogram %s",runNumbers[ri],Form("run%i_hPi0MassSM%iSM%i",runNumbers[ri],SM1,SM1)); | |
1931 | continue; | |
1932 | } | |
9f3c053c | 1933 | |
1934 | // correct for SM acceptance | |
1935 | if (hmap) { | |
1936 | Double_t accep = GetAcceptance(sm1, hmap, ri); | |
1937 | if (accep > 0) h->Scale(1./accep); | |
1938 | else noSM++; | |
1939 | } | |
1940 | ||
1941 | hsum->Add(h); | |
1942 | delete h; | |
1943 | } | |
1944 | ||
1945 | Double_t nraw, enraw, mass, emass, sigma, esigma; | |
1946 | FitPi0(hsum, nraw, enraw, mass, emass, sigma, esigma); | |
1947 | ||
35ac8df7 | 1948 | hPi0Num ->SetBinContent(ri+1, nraw/nevents/(SM2-SM1+1-noSM)); |
1949 | hPi0Mass ->SetBinContent(ri+1, mass); | |
9f3c053c | 1950 | hPi0Sigma->SetBinContent(ri+1, sigma); |
1951 | ||
35ac8df7 | 1952 | hPi0Num ->SetBinError(ri+1, enraw/nevents/(SM2-SM1+1-noSM)); |
1953 | hPi0Mass ->SetBinError(ri+1, emass); | |
9f3c053c | 1954 | hPi0Sigma->SetBinError(ri+1, esigma); |
1955 | ||
1956 | delete hsum; | |
1957 | } // run index loop | |
1958 | ||
1959 | fprintf(stderr, "\n"); | |
1960 | ||
1961 | ||
1962 | /* Draw results | |
1963 | */ | |
1964 | ||
1965 | // number of pi0s vs run index | |
1966 | TCanvas *c1 = new TCanvas(Form("hPi0NumRuns%s",suff),Form("hPi0NumRuns%s",suff), 800,400); | |
1967 | gPad->SetLeftMargin(0.08); | |
c9365469 | 1968 | gPad->SetRightMargin(0.03); |
9f3c053c | 1969 | gPad->SetTopMargin(0.12); |
35ac8df7 | 1970 | gPad->SetBottomMargin(0.12); |
1971 | gPad->SetGridx(); | |
1972 | gPad->SetGridy(); | |
1973 | // TLegend *leg = new TLegend(0.75,0.15,0.89,0.15+0.06*(SM2-SM1+1)); | |
1974 | TLegend *leg = new TLegend(0.75,0.89-0.06*(SM2-SM1+1),0.89,0.89); | |
9f3c053c | 1975 | |
35ac8df7 | 1976 | hPi0Num->SetAxisRange(0., hPi0Num->GetMaximum()*1.8, "Y"); |
c9365469 | 1977 | hPi0Num->SetTitleOffset(0.8, "Y"); |
1978 | hPi0Num->SetLineWidth(2); | |
35ac8df7 | 1979 | hPi0Num->SetLabelSize(0.04,"X"); |
9f3c053c | 1980 | hPi0Num->Draw(); |
c9365469 | 1981 | leg->AddEntry(hPi0Num, Form("(All Modules)/%i",SM2-SM1+1),"l"); |
9f3c053c | 1982 | for (Int_t sm = SM1; sm <= SM2; sm++) { |
1983 | hPi0NumSM[sm]->SetLineColor(colorsSM[sm]); | |
c9365469 | 1984 | hPi0NumSM[sm]->SetLineWidth(2); |
9f3c053c | 1985 | hPi0NumSM[sm]->Draw("same"); |
c9365469 | 1986 | leg->AddEntry(hPi0NumSM[sm], Form("Module %i",5-sm),"l"); |
9f3c053c | 1987 | } |
1988 | hPi0Num->Draw("same"); // to the top | |
1989 | hPi0Num->Draw("hist,same"); | |
1990 | leg->Draw("same"); | |
35ac8df7 | 1991 | c1->Update(); |
9f3c053c | 1992 | |
35ac8df7 | 1993 | if (trendPlotEvents) { |
1994 | ||
1995 | // number of pi0s vs event count | |
1996 | TCanvas *c2 = new TCanvas(Form("hPi0NumEvents%s",suff),Form("hPi0NumEvents%s",suff), 800,400); | |
1997 | gPad->SetLeftMargin(0.08); | |
1998 | gPad->SetRightMargin(0.06); | |
1999 | gPad->SetTopMargin(0.12); | |
2000 | leg = new TLegend(0.75,0.15,0.92,0.15+0.04*(SM2-SM1+1)); | |
2001 | ||
2002 | TH1* hPi0NumEvents = ConvertMapRuns2Events(nruns, runNumbers, hPi0Num); | |
2003 | hPi0NumEvents->Draw(); | |
2004 | leg->AddEntry(hPi0NumEvents, Form("(All Modules)/%i",SM2-SM1+1),"l"); | |
2005 | for (Int_t sm = SM1; sm <= SM2; sm++) { | |
2006 | TH1* hPi0NumSMEvents = ConvertMapRuns2Events(nruns, runNumbers, hPi0NumSM[sm]); | |
2007 | hPi0NumSMEvents->Draw("same"); | |
2008 | leg->AddEntry(hPi0NumSMEvents, Form("Module %i",5-sm),"l"); | |
2009 | } | |
2010 | hPi0NumEvents->Draw("same"); // to the top | |
2011 | hPi0NumEvents->Draw("hist,same"); | |
2012 | leg->Draw("same"); | |
2013 | c2->Update(); | |
9f3c053c | 2014 | } |
35ac8df7 | 2015 | |
9f3c053c | 2016 | // pi0 mass and width vs run index |
c9365469 | 2017 | TCanvas *c3 = new TCanvas(Form("hPi0MassSigmaRuns%s",suff),Form("hPi0MassSigmaRuns%s",suff), 1000,707); |
2018 | c3->Divide(1,2); | |
9f3c053c | 2019 | |
2020 | c3->cd(1); | |
c9365469 | 2021 | gPad->SetLeftMargin(0.06); |
9f3c053c | 2022 | gPad->SetRightMargin(0.04); |
c9365469 | 2023 | gPad->SetTopMargin(0.10); |
2024 | gPad->SetBottomMargin(0.13); | |
35ac8df7 | 2025 | gPad->SetGridx(); |
2026 | gPad->SetGridy(); | |
2027 | // leg = new TLegend(0.75,0.15,0.85,0.16+0.06*(SM2-SM1+1)); | |
2028 | leg = new TLegend(0.75,0.89-0.06*(SM2-SM1+1),0.85,0.89); | |
9f3c053c | 2029 | |
35ac8df7 | 2030 | hPi0Mass->SetAxisRange(0.125, 0.149, "Y"); |
2031 | hPi0Mass->SetTitleOffset(0.8, "Y"); | |
c9365469 | 2032 | hPi0Mass->SetLineWidth(2); |
35ac8df7 | 2033 | hPi0Mass->SetLabelSize(0.06,"X"); |
9f3c053c | 2034 | hPi0Mass->Draw(); |
c9365469 | 2035 | leg->AddEntry(hPi0Mass, "All Modules","l"); |
9f3c053c | 2036 | for (Int_t sm = SM1; sm <= SM2; sm++) { |
2037 | hPi0MassSM[sm]->SetLineColor(colorsSM[sm]); | |
c9365469 | 2038 | hPi0MassSM[sm]->SetLineWidth(2); |
9f3c053c | 2039 | hPi0MassSM[sm]->Draw("same"); |
c9365469 | 2040 | leg->AddEntry(hPi0MassSM[sm], Form("Module %i",5-sm),"l"); |
9f3c053c | 2041 | } |
2042 | hPi0Mass->Draw("same"); // to the top | |
2043 | hPi0Mass->Draw("hist,same"); | |
2044 | leg->Draw("same"); | |
2045 | ||
2046 | c3->cd(2); | |
c9365469 | 2047 | gPad->SetLeftMargin(0.06); |
9f3c053c | 2048 | gPad->SetRightMargin(0.04); |
c9365469 | 2049 | gPad->SetTopMargin(0.10); |
2050 | gPad->SetBottomMargin(0.13); | |
35ac8df7 | 2051 | gPad->SetGridx(); |
2052 | gPad->SetGridy(); | |
2053 | // leg = new TLegend(0.75,0.15,0.85,0.16+0.06*(SM2-SM1+1)); | |
2054 | leg = new TLegend(0.75,0.89-0.06*(SM2-SM1+1),0.85,0.89); | |
9f3c053c | 2055 | |
2056 | hPi0Sigma->SetAxisRange(0., hPi0Sigma->GetMaximum()*1.5, "Y"); | |
c9365469 | 2057 | hPi0Sigma->SetTitleOffset(0.8, "Y"); |
2058 | hPi0Sigma->SetLineWidth(2); | |
35ac8df7 | 2059 | hPi0Sigma->SetLabelSize(0.06,"X"); |
9f3c053c | 2060 | hPi0Sigma->Draw(); |
c9365469 | 2061 | leg->AddEntry(hPi0Sigma, "All Modules","l"); |
9f3c053c | 2062 | for (Int_t sm = SM1; sm <= SM2; sm++) { |
2063 | hPi0SigmaSM[sm]->SetLineColor(colorsSM[sm]); | |
c9365469 | 2064 | hPi0SigmaSM[sm]->SetLineWidth(2); |
9f3c053c | 2065 | hPi0SigmaSM[sm]->Draw("same"); |
c9365469 | 2066 | leg->AddEntry(hPi0SigmaSM[sm], Form("Module %i",5-sm),"l"); |
9f3c053c | 2067 | } |
2068 | hPi0Sigma->Draw("same"); // to the top | |
2069 | hPi0Sigma->Draw("hist,same"); | |
2070 | leg->Draw("same"); | |
35ac8df7 | 2071 | c3->Update(); |
9f3c053c | 2072 | |
2073 | ||
35ac8df7 | 2074 | if (trendPlotEvents) { |
2075 | // pi0 mass and width vs number of events | |
2076 | TCanvas *c4 = new TCanvas(Form("hPi0MassSigmaEvents%s",suff),Form("hPi0MassSigmaEvents%s",suff), 800,400); | |
2077 | c4->Divide(1,2); | |
2078 | ||
2079 | c4->cd(1); | |
2080 | gPad->SetLeftMargin(0.16); | |
2081 | gPad->SetRightMargin(0.08); | |
2082 | leg = new TLegend(0.75,0.15,0.91,0.15+0.04*(SM2-SM1+1)); | |
2083 | ||
2084 | TH1* hPi0MassEvents = ConvertMapRuns2Events(nruns, runNumbers, hPi0Mass); | |
2085 | hPi0MassEvents->Draw(); | |
2086 | leg->AddEntry(hPi0MassEvents, "All Modules","l"); | |
2087 | for (Int_t sm = SM1; sm <= SM2; sm++) { | |
2088 | TH1* hPi0MassSMEvents = ConvertMapRuns2Events(nruns, runNumbers, hPi0MassSM[sm]); | |
2089 | hPi0MassSMEvents->Draw("same"); | |
2090 | leg->AddEntry(hPi0MassSMEvents, Form("Module %i",5-sm),"l"); | |
2091 | } | |
2092 | hPi0MassEvents->Draw("same"); // to the top | |
2093 | hPi0MassEvents->Draw("hist,same"); | |
2094 | leg->Draw("same"); | |
2095 | ||
2096 | c4->cd(2); | |
2097 | gPad->SetLeftMargin(0.16); | |
2098 | gPad->SetRightMargin(0.08); | |
2099 | leg = new TLegend(0.75,0.15,0.91,0.15+0.04*(SM2-SM1+1)); | |
2100 | ||
2101 | TH1* hPi0SigmaEvents = ConvertMapRuns2Events(nruns, runNumbers, hPi0Sigma); | |
2102 | hPi0SigmaEvents->Draw(); | |
2103 | leg->AddEntry(hPi0SigmaEvents, "All Modules","l"); | |
2104 | for (Int_t sm = SM1; sm <= SM2; sm++) { | |
2105 | TH1* hPi0SigmaSMEvents = ConvertMapRuns2Events(nruns, runNumbers, hPi0SigmaSM[sm]); | |
2106 | hPi0SigmaSMEvents->Draw("same"); | |
2107 | leg->AddEntry(hPi0SigmaSMEvents, Form("Module %i",5-sm),"l"); | |
2108 | } | |
2109 | hPi0SigmaEvents->Draw("same"); // to the top | |
2110 | hPi0SigmaEvents->Draw("hist,same"); | |
2111 | leg->Draw("same"); | |
2112 | c4->Update(); | |
9f3c053c | 2113 | } |
9f3c053c | 2114 | |
9f3c053c | 2115 | } |
2116 | ||
2117 | //_________________________________________________________________________ | |
2118 | void DrawPi0Distributions(char *suff, Int_t nbins = 100) | |
2119 | { | |
2120 | // Draw distributions for | |
2121 | // 1) average number of pi0s per event; | |
2122 | // 2) pi0 mass position; | |
2123 | // 3) pi0 width. | |
2124 | // | |
2125 | // Must be called after DrawPi0Averages() because it | |
2126 | // searches for the corresponding histograms by name. | |
2127 | // | |
2128 | // suff -- histograms suffix; | |
2129 | // nbins -- number of bins in distributions. | |
2130 | ||
2131 | TCanvas *c1 = new TCanvas(Form("Pi0Distributions%s",suff),Form("Pi0Distributions%s",suff), 999,333); | |
2132 | c1->Divide(3,1); | |
2133 | ||
2134 | // number of pi0s | |
2135 | c1->cd(1)->SetLogy(); | |
2136 | gPad->SetLeftMargin(0.16); | |
2137 | gPad->SetRightMargin(0.04); | |
2138 | TH1* hPi0Num = (TH1*) gROOT->FindObject(Form("hPi0Num%s",suff)); | |
2139 | MakeDistribution(hPi0Num,nbins)->Draw(); | |
2140 | ||
2141 | // pi0 mass | |
2142 | c1->cd(2)->SetLogy(); | |
2143 | gPad->SetLeftMargin(0.16); | |
2144 | gPad->SetRightMargin(0.04); | |
2145 | TH1* hPi0Mass = (TH1*) gROOT->FindObject(Form("hPi0Mass%s",suff)); | |
2146 | MakeDistribution(hPi0Mass,nbins)->Draw(); | |
2147 | ||
2148 | // pi0 width | |
2149 | c1->cd(3)->SetLogy(); | |
2150 | gPad->SetLeftMargin(0.16); | |
2151 | gPad->SetRightMargin(0.04); | |
2152 | TH1* hPi0Sigma = (TH1*) gROOT->FindObject(Form("hPi0Sigma%s",suff)); | |
2153 | MakeDistribution(hPi0Sigma,nbins)->Draw(); | |
2154 | ||
2155 | c1->Update(); | |
2156 | } | |
2157 | ||
2158 | //_________________________________________________________________________ | |
2159 | TH1* MakeDistribution(TH1* inhisto, Int_t nbins) | |
2160 | { | |
2161 | // Returns distribution for the input histogram. | |
2162 | // | |
2163 | // inhisto -- input histogram; | |
2164 | // nbins -- number of bins in distribution. | |
2165 | ||
2166 | // initialize distribution; 1.01 -- to see the last bin | |
2167 | TH1* distr = new TH1F(Form("%sDistr",inhisto->GetName()),"", | |
2168 | nbins, inhisto->GetMinimum(),inhisto->GetMaximum()*1.01); | |
2169 | distr->SetXTitle(inhisto->GetYaxis()->GetTitle()); | |
2170 | distr->SetYTitle("Number of runs"); | |
2171 | ||
2172 | // fill distribution | |
2173 | for (Int_t i = 1; i <= inhisto->GetNbinsX(); i++) | |
2174 | distr->Fill(inhisto->GetBinContent(i)); | |
2175 | ||
2176 | return distr; | |
2177 | } | |
2178 | ||
2179 | //_________________________________________________________________________ | |
2180 | void FitPi0(TH1* h, Double_t &nraw, Double_t &enraw, | |
2181 | Double_t &mass, Double_t &emass, | |
2182 | Double_t &sigma, Double_t &esigma, | |
2183 | Double_t emin = 0.05, Double_t emax = 0.3, Int_t rebin = 1) | |
2184 | { | |
2185 | // Fits the pi0 peak with crystal ball + pol2, | |
2186 | // fills number of pi0s, mass, width and their errors. | |
2187 | ||
2188 | nraw = enraw = 0; | |
2189 | mass = emass = 0; | |
2190 | sigma = esigma = 0; | |
2191 | ||
2192 | if (h->GetEntries() == 0) return NULL; | |
2193 | ||
2194 | if (rebin > 1) h->Rebin(rebin); | |
2195 | ||
2196 | // crystal ball parameters (fixed by hand for EMCAL); TODO: PHOS parameters? | |
2197 | Double_t alpha = 1.1; // alpha >= 0 | |
2198 | Double_t n = 2.; // n > 1 | |
2199 | ||
2200 | // CB tail parameters | |
2201 | Double_t a = pow((n/alpha), n) * TMath::Exp(-alpha*alpha/2.); | |
2202 | Double_t b = n/alpha - alpha; | |
2203 | ||
2204 | // integral value under crystal ball with amplitude = 1, sigma = 1 | |
2205 | // (will be sqrt(2pi) at alpha = infinity) | |
2206 | Double_t nraw11 = a * pow(b+alpha, 1.-n)/(n-1.) + TMath::Sqrt(TMath::Pi()/2.) * TMath::Erfc(-alpha/TMath::Sqrt(2.)); | |
2207 | ||
2208 | // signal (crystal ball); | |
2209 | new TF1("cball", Form("(x-[1])/[2] > -%f ? \ | |
2210 | [0]*exp(-(x-[1])*(x-[1])/(2*[2]*[2])) \ | |
2211 | : [0]*%f*(%f-(x-[1])/[2])^(-%f)", alpha, a, b, n)); | |
2212 | ||
2213 | // background | |
2214 | new TF1("mypol2", "[0] + [1]*(x-0.135) + [2]*(x-0.135)^2", emin, emax); | |
2215 | ||
2216 | // signal + background | |
2217 | TF1 *fitfun = new TF1("fitfun", "cball + mypol2", emin, emax); | |
2218 | fitfun->SetParNames("A", "M", "#sigma", "a_{0}", "a_{1}", "a_{2}"); | |
2219 | fitfun->SetLineColor(kRed); | |
2220 | fitfun->SetLineWidth(2); | |
2221 | ||
2222 | // make a preliminary fit to estimate parameters | |
2223 | TF1* ff = new TF1("fastfit", "gaus(0) + [3]"); | |
35ac8df7 | 2224 | ff->SetParLimits(0, 0., h->GetMaximum()*1.5); |
9f3c053c | 2225 | ff->SetParLimits(1, 0.1, 0.2); |
35ac8df7 | 2226 | ff->SetParLimits(2, 0.004,0.030); |
9f3c053c | 2227 | ff->SetParameters(h->GetMaximum()/3., 0.135, 0.010, 0.); |
35ac8df7 | 2228 | h->Fit(ff, "0QL", "", 0.105, 0.165); |
9f3c053c | 2229 | |
35ac8df7 | 2230 | fitfun->SetParLimits(0, 0., h->GetMaximum()*1.5); |
2231 | fitfun->SetParLimits(1, 0.12, 0.15); | |
2232 | fitfun->SetParLimits(2, 0.004,0.030); | |
9f3c053c | 2233 | fitfun->SetParameters(ff->GetParameter(0), ff->GetParameter(1), ff->GetParameter(2), ff->GetParameter(3)); |
35ac8df7 | 2234 | h->Fit(fitfun,"QLR", ""); |
9f3c053c | 2235 | |
2236 | // calculate pi0 values | |
2237 | mass = fitfun->GetParameter(1); | |
2238 | emass = fitfun->GetParError(1); | |
2239 | ||
2240 | sigma = fitfun->GetParameter(2); | |
2241 | esigma = fitfun->GetParError(2); | |
2242 | ||
2243 | Double_t A = fitfun->GetParameter(0); | |
2244 | Double_t eA = fitfun->GetParError(0); | |
2245 | ||
2246 | nraw = nraw11 * A * sigma / h->GetBinWidth(1); | |
2247 | enraw = nraw * (eA/A + esigma/sigma); | |
2248 | } | |
2249 | ||
2250 | ||
2251 | /* COMMON FUNCTIONS FOR RUN NUMBERS EXTRACTION, VISUALIZATION AND FILTERING. | |
2252 | * | |
2253 | * Input: histogram hNEventsProcessedPerRun which is taken from gInputFile. | |
2254 | */ | |
2255 | ||
2256 | //_________________________________________________________________________ | |
2257 | void GetRunNumbers(Int_t &nruns, Int_t runNumbers[]) | |
2258 | { | |
2259 | // Fill runNumbers array with run numbers according to hNEventsProcessedPerRun | |
2260 | // histogram: actually analysed runs have positive bin content. | |
2261 | ||
2262 | TH1* hNEventsProcessedPerRun = (TH1*) gInputFile->Get("hNEventsProcessedPerRun"); | |
2263 | ||
2264 | // check underflow/overflow | |
2265 | if (hNEventsProcessedPerRun->GetBinContent(0) > 0.5) | |
2266 | Warning("GetRunNumbers", "some run numbers are in underflow; they will be absent in the list of runs"); | |
2267 | if (hNEventsProcessedPerRun->GetBinContent(hNEventsProcessedPerRun->GetNbinsX()+1) > 0.5) | |
2268 | Warning("GetRunNumbers", "some run numbers are in overflow; they will be absent in the list of runs"); | |
2269 | ||
2270 | nruns = 0; | |
2271 | ||
2272 | for (Int_t i = 1; i <= hNEventsProcessedPerRun->GetNbinsX(); i++) | |
2273 | if (hNEventsProcessedPerRun->GetBinContent(i) > 0.5) { | |
2274 | runNumbers[nruns] = hNEventsProcessedPerRun->GetBinLowEdge(i); | |
2275 | nruns++; | |
2276 | } | |
2277 | } | |
2278 | ||
2279 | //_________________________________________________________________________ | |
2280 | Int_t GetNumberOfEvents(Int_t run) | |
2281 | { | |
2282 | // Return number of events in run; | |
2283 | // run -- run number. | |
2284 | ||
2285 | TH1* hNEventsProcessedPerRun = (TH1*) gInputFile->Get("hNEventsProcessedPerRun"); | |
2286 | ||
2287 | // round the number of events to avoid precision surprizes | |
2288 | return TMath::Nint( hNEventsProcessedPerRun->GetBinContent(hNEventsProcessedPerRun->FindBin(run)) ); | |
2289 | } | |
2290 | ||
2291 | //_________________________________________________________________________ | |
2292 | Long64_t GetTotalNumberOfEvents(Int_t nruns, Int_t runNumbers[]) | |
2293 | { | |
2294 | // Return total number of events in all runs | |
2295 | ||
2296 | Long64_t ntotal = 0; | |
2297 | ||
2298 | for (Int_t i = 0; i < nruns; i++) | |
2299 | ntotal += GetNumberOfEvents(runNumbers[i]); | |
2300 | ||
2301 | return ntotal; | |
2302 | } | |
2303 | ||
2304 | //_________________________________________________________________________ | |
2305 | void DrawRunsDistribution(Int_t nruns, Int_t runNumbers[], Int_t dnbins = 100) | |
2306 | { | |
2307 | // Draw events and events distribution; | |
2308 | // dnbins -- number of bins in distribution. | |
2309 | ||
2310 | // initialize run index histogram ... | |
2311 | TH1* hNEventsPerRunIndex = new TH1F("hNEventsPerRunIndex", "Number of processed events per run index", nruns,0,nruns); | |
c9365469 | 2312 | SetRunLabel(hNEventsPerRunIndex,nruns,runNumbers,1); |
2313 | // hNEventsPerRunIndex->SetXTitle("Run index"); | |
9f3c053c | 2314 | hNEventsPerRunIndex->SetYTitle("Number of events"); |
2315 | ||
2316 | // ... and fill it | |
2317 | for (Int_t i = 0; i < nruns; i++) | |
2318 | hNEventsPerRunIndex->SetBinContent(i+1, GetNumberOfEvents(runNumbers[i])); | |
2319 | ||
2320 | // initialize distribution; 1.01 - to see the last bin | |
2321 | TH1* distrib = new TH1F("hNEventsPerRunIndexDistr", "", dnbins, 0, hNEventsPerRunIndex->GetMaximum()*1.01); | |
2322 | distrib->SetXTitle("Number of events"); | |
2323 | distrib->SetYTitle("Number of runs"); | |
2324 | ||
2325 | // fill distribution | |
2326 | for (Int_t i = 1; i <= nruns; i++) | |
2327 | distrib->Fill(hNEventsPerRunIndex->GetBinContent(i)); | |
2328 | ||
2329 | // draw histogram + distribution | |
c9365469 | 2330 | TCanvas *c1 = new TCanvas("hNEventsPerRunIndex","hNEventsPerRunIndex", 800,600); |
2331 | c1->Divide(1,2); | |
9f3c053c | 2332 | |
2333 | c1->cd(1); | |
c9365469 | 2334 | gPad->SetLeftMargin(0.06); |
2335 | gPad->SetRightMargin(0.04); | |
2336 | gPad->SetTopMargin(0.10); | |
35ac8df7 | 2337 | gPad->SetBottomMargin(0.14); |
2338 | gPad->SetGridx(); | |
c9365469 | 2339 | hNEventsPerRunIndex->SetTitleOffset(0.6,"Y"); |
2340 | hNEventsPerRunIndex->SetTickLength(0.01,"Y"); | |
35ac8df7 | 2341 | hNEventsPerRunIndex->SetLabelSize(0.06,"X"); |
9f3c053c | 2342 | hNEventsPerRunIndex->Draw(); |
2343 | ||
2344 | c1->cd(2); | |
c9365469 | 2345 | gPad->SetLeftMargin(0.06); |
2346 | gPad->SetRightMargin(0.04); | |
2347 | gPad->SetTopMargin(0.10); | |
2348 | gPad->SetBottomMargin(0.13); | |
2349 | distrib->SetTitleOffset(0.6,"Y"); | |
2350 | distrib->SetTickLength(0.01,"Y"); | |
9f3c053c | 2351 | distrib->Draw(); |
2352 | ||
2353 | c1->Update(); | |
2354 | } | |
2355 | ||
2356 | //_________________________________________________________________________ | |
2357 | void ExcludeSmallRuns(Int_t &nruns, Int_t runNumbers[], Int_t nEventsMin = 100e+3) | |
2358 | { | |
2359 | // Exclude runs with number of events < nEventsMin | |
2360 | ||
2361 | printf("Excluding runs (< %.0fk events):", nEventsMin/1000.); | |
2362 | ||
2363 | for (Int_t i = 0; i < nruns; i++) | |
2364 | if (GetNumberOfEvents(runNumbers[i]) < nEventsMin) { | |
2365 | printf(" %i", runNumbers[i]); | |
2366 | nruns--; | |
2367 | runNumbers[i] = runNumbers[nruns]; | |
2368 | i--; | |
2369 | } | |
2370 | ||
2371 | printf("\n"); | |
2372 | ||
2373 | SortArray(nruns, runNumbers); | |
2374 | } | |
2375 | ||
2376 | //_________________________________________________________________________ | |
2377 | void ExcludeRunNumbers(Int_t &nruns, Int_t runNumbers[], Int_t nexcl, Int_t runs2Exclude[]) | |
2378 | { | |
2379 | // Exclude particular runs. | |
2380 | // | |
2381 | // nexcl -- number of runs in runs2Exclude[]; | |
2382 | // runs2Exclude -- array with run numbers to exclude. | |
2383 | ||
2384 | for (Int_t i = 0; i < nruns; i++) | |
2385 | for (Int_t e = 0; e < nexcl; e++) | |
2386 | if (runNumbers[i] == runs2Exclude[e]) { | |
2387 | Printf("Excluding run: %i", runs2Exclude[e]); | |
2388 | nruns--; | |
2389 | runNumbers[i] = runNumbers[nruns]; | |
2390 | i--; | |
2391 | break; | |
2392 | } | |
2393 | ||
2394 | SortArray(nruns, runNumbers); | |
2395 | } | |
2396 | ||
2397 | //_________________________________________________________________________ | |
2398 | void SortArray(const Int_t nruns, Int_t runNumbers[]) | |
2399 | { | |
2400 | // Sort an array; | |
2401 | // used for runNumbers array sorting. | |
2402 | ||
2403 | Int_t indices[nruns]; | |
2404 | Int_t runNumbers_unsort[nruns]; | |
2405 | ||
2406 | for (Int_t i = 0; i < nruns; i++) | |
2407 | runNumbers_unsort[i] = runNumbers[i]; | |
2408 | ||
2409 | TMath::Sort(nruns, runNumbers_unsort, indices, kFALSE); | |
2410 | ||
2411 | for (Int_t i = 0; i < nruns; i++) | |
2412 | runNumbers[i] = runNumbers_unsort[indices[i]]; | |
2413 | } | |
2414 | ||
c9365469 | 2415 | //_________________________________________________________________________ |
9f3c053c | 2416 | void PrintRunNumbers(Int_t nruns, Int_t runNumbers[]) |
2417 | { | |
2418 | // Print a table run index / run number / number of events. | |
2419 | ||
2420 | Printf(""); | |
2421 | Printf("| index | run number | nevents |"); | |
2422 | Printf("|-------|------------|-----------|"); | |
2423 | ||
2424 | for (Int_t i = 0; i < nruns; i++) | |
2425 | Printf("| %-4i | %-9i | %-8i |", i, runNumbers[i], GetNumberOfEvents(runNumbers[i])); | |
2426 | ||
2427 | Printf("| Events in total: %-10lli |\n", GetTotalNumberOfEvents(nruns, runNumbers)); | |
2428 | } | |
c9365469 | 2429 | |
2430 | //------------------------------------------------------------------------ | |
2431 | void SetRunLabel(TH1 *histo, Int_t nruns, Int_t *runNumbers, Int_t axis) | |
2432 | { | |
2433 | TString runText; | |
2434 | for (Int_t i=0; i<nruns; i++) { | |
2435 | runText = Form("%d",runNumbers[i]); | |
2436 | if (axis == 1) | |
2437 | histo->GetXaxis()->SetBinLabel(i+1,runText.Data()); | |
2438 | else if (axis == 2) | |
2439 | histo->GetYaxis()->SetBinLabel(i+1,runText.Data()); | |
2440 | } | |
2441 | histo->LabelsOption("v"); | |
2442 | } |