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