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.
9 * For impatient, change the "infile" line below and execute
10 * root -l getCellsRunsQA.C
12 * For curios, continue reading getCellsRunsQA() code: it is self-documented.
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.
25 * Generally, a QA expert uncomments all the functions which return (print to stdout)
26 * bad cell candidates and checks them by hand.
28 * Detector-specific parts require to run this macro with aliroot instead of root,
29 * they are commented in getCellsRunsQA() by default.
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.
35 * Input: AliCaloCellsQA analysis output.
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...
41 * TODO: some PHOS-specific parts
43 * Author: Olga Driga (SUBATECH)
47 // ========================== Parameters section ==============================
50 char *infile = "CellsQAEMCAL.root";
53 Color_t colorsSM[] = {2,3,4,5,6,7,8,9,11,12};
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 Int_t nexc = 45+69+156+132;
60 74,103,152,495,871,917,1059,1263,1275,1276,1288,1376,1382,1384,1519,1712,1961, // EMCAL, LHC10bcde
61 1967,2026,2047,2112,2114,2115,2116,2117,2120,2123,2298,2540,2671,2768,2769,2770, // 45
62 2771,2773,2774,2776,2777,2778,2779,2780,2783,3135,3544,3567
64 ,49,152,495,917,1248,1249,1250,1251,1252,1253,1254,1255,1256,1257,1258,1259,1260, // EMCAL, LHC10h
65 1261,1262,1263,1275,1276,1288,1296,1297,1298,1299,1300,1301,1302,1303,1304,1305, // 69
66 1306,1307,1308,1309,1310,1311,1376,1384,1414,1519,1712,1860,1967,2026,2047,2114,
67 2115,2116,2117,2120,2123,2298,2671,2768,2769,2770,2771,2773,2774,2776,2777,2778,
70 ,74,103,320,321,322,323,324,325,326,327,328,329,330,331,332,333,334,335,368,369, // EMCAL, LHC11a pp@7TeV
71 370,371,372,373,374,375,376,377,378,379,380,381,382,383,1263,1275,1276,1288,1376, // 156
72 1382,1384,1519,1595,1704,1860,1967,2022,2026,2047,2071,2114,2115,2116,2117,2120,
73 2123,2210,2298,2776,2778,3544,3567,4026,4157,4174,4494,4609,4610,4616,4620,4662,
74 4670,4705,4718,4719,4752,4754,4756,4758,4762,4764,4800,4801,4802,4804,4857,4898,
75 4904,4958,4997,4999,5041,5043,5044,5090,5097,5102,5137,5138,5139,5148,5184,5192,
76 5208,5246,5330,5336,5384,5428,5430,5436,5437,5726,5748,5754,5767,6095,6111,6340,
77 6592,6800,6801,6802,6803,6804,6805,6806,6807,6808,6809,6810,6811,6812,6813,6814,
78 6815,7089,7371,7425,7430,7457,7491,7709,7874,8352,8353,8356,8357,8808,8810,8812,
79 8814,9056,9815,9837,11080
81 ,74,103,320,321,322,323,324,325,326,327,328,329,330,331,332,333,334,335,368,369, // EMCAL, LHC11a pp@2.76TeV
82 370,371,372,373,374,375,376,377,378,379,380,381,382,383,1059,1175,1263,1275,1276, // 132
83 1288,1376,1382,1384,1519,1595,1860,1967,2022,2026,2047,2071,2117,2210,2298,2362,
84 2506,2776,2778,3544,4026,4157,4174,4494,4610,4616,4620,4662,4670,4719,4756,4758,
85 4764,4801,4857,4898,4904,4958,4999,5041,5090,5097,5102,5138,5148,5184,5192,5246,
86 5330,5336,5384,5430,5437,5767,6095,6111,6340,6592,6800,6801,6802,6803,6804,6805,
87 6806,6807,6808,6809,6810,6811,6812,6813,6814,6815,7371,7375,7425,7430,7457,7491,
88 7709,7874,8352,8353,8356,8357,8808,8810,8812,8814,9056,9815,9837
91 // ====================== End of parameters section ===========================
96 // Entry point for the analysis.
99 gStyle->SetOptStat(0);
100 gStyle->SetOptFit(111);
101 gStyle->SetPalette(1);
102 gStyle->SetFillColor(10);
104 Printf("Input: %s", infile);
107 // You may wish to extract and draw a particular histogram from infile.
108 // Here are the examples:
110 // TH1* hNEventsProcessedPerRun = (TH1*) gFile->Get("hNEventsProcessedPerRun");
111 // hNEventsProcessedPerRun->Draw();
114 // TH1* hNTimesInClusterElow = (TH1*) gFile->Get("run146686_hCellLocMaxNTimesInClusterElow");
115 // hNTimesInClusterElow->Draw();
118 // TH1* hETotalClusterEhigh = (TH1*) gFile->Get("run146686_hCellLocMaxETotalClusterEhigh");
119 // hETotalClusterEhigh->Draw();
122 // Draw a random cell spectrum;
123 // 0.2-1 GeV -- fit region;
124 // 3 GeV -- the spectrum is drawn up to this value, -1 = no limit;
125 // last argument -- histogram name to process, possible values are:
126 // hCellAmplitude, hCellAmplitudeEHigh, hCellAmplitudeNonLocMax or fhCellAmplitudeEhighNonLocMax.
127 DrawCell(1+gRandom->Integer(4607), 0.25, 1., 2., "hCellAmplitude"); // common cell region for EMCAL2010 and for PHOS
129 // Draw a random cell time spectrum
130 DrawCellTime(1+gRandom->Integer(4607));
133 /* RUN NUMBER SELECTION SECTION
135 * NOTE: at any point below runs are sorted in chronological order.
138 // array with run numbers and their number
139 Int_t runNumbers[10000];
142 // First, fill run numbers ...
143 GetRunNumbers(nruns, runNumbers);
144 Printf("Total number of runs: %i", nruns);
146 // ... draw events distribution ...
147 // (the last argument is number of bins in this distribution)
148 DrawRunsDistribution(nruns, runNumbers, 100);
150 // ... and exclude runs with number of events < 50k.
151 ExcludeSmallRuns(nruns, runNumbers, 100e+3);
153 // You may wish to exclude particular runs:
154 // Int_t runs2Exclude[] = {111222,333444,555666};
155 // ExcludeRunNumbers(nruns, runNumbers, 3, runs2Exclude);
157 Printf("Number of runs to be analysed: %i", nruns);
159 // Finally, print a nice table with run index / run number / number of events.
160 PrintRunNumbers(nruns, runNumbers);
163 /* PER RUN BAD CELLS SEARCHING CRITERIA
165 * Four primary criteria on a per run basis:
166 * 1 and 2: number of times cell was a local maximum in a cluster at low/high energies;
167 * 3 and 4: total cluster energy for a local maximum cell at low/high energies.
170 // Extract the histograms with dead/noisy cell candidates per run.
171 // For each of the four criteria, for each run and for each cell:
172 // 1) calculate cell factor = [cell value]/[average over cells];
173 // 2) mark cell as dead (candidate) if factor <= factorDead (3rd argument below);
174 // 3) mark cell as noisy (candidate) if factor >= factorNoisy (4th argument below).
176 // Factor thresholds are quite wide by default:
177 // low energy criteria are not very sensitive to them,
178 // while high energy criteria are very sensitive due to limited statistics.
180 // The function below also draws histograms with factor distributions in all the runs
181 // for all the cells. It may help to take the decision about dead/noisy factor thresholds.
182 // The last two arguments -- number of bins and maximum X axis value for such histograms.
184 // Output: hBadCellMap[4] -- bad cell candidates per run for each of the four criteria;
185 // X axis -- cell absId, Y axis -- run index, content: 0=not bad, 1=noisy, -1=dead.
187 TH2** hBadCellMap = FindDeadNoisyCellsPerRun(nruns, runNumbers, 0.05, 3.5, 200, 10.);
189 // Print cell numbers suggested for exclusion from averages calculation;
190 // see excells[] array in the parameters section in the beginning of the macro.
191 // It is highly suggested to run this macro several times and add the output
192 // of this function to the parameters section.
193 PrintCellsToExcludeFromAverages(hBadCellMap);
195 // The results can be visualized (second argument -- canvas name):
196 // either per run ...
197 DrawDeadNoisyCellMap(hBadCellMap[0], "hMapRuns_NTimesInClusterElow");
198 DrawDeadNoisyCellMap(hBadCellMap[1], "hMapRuns_NTimesInClusterEhigh");
199 // DrawDeadNoisyCellMap(hBadCellMap[2], "hMapRuns_ETotalClusterElow");
200 // DrawDeadNoisyCellMap(hBadCellMap[3], "hMapRuns_ETotalClusterEhigh");
202 // ... or per number of events ...
203 // DrawDeadNoisyCellMap((TH2*)ConvertMapRuns2Events(nruns,runNumbers,hBadCellMap[0]), "hMapEvents_NTimesInClusterElow");
204 // DrawDeadNoisyCellMap((TH2*)ConvertMapRuns2Events(nruns,runNumbers,hBadCellMap[1]), "hMapEvents_NTimesInClusterEhigh");
205 // DrawDeadNoisyCellMap((TH2*)ConvertMapRuns2Events(nruns,runNumbers,hBadCellMap[2]), "hMapEvents_ETotalClusterElow");
206 // DrawDeadNoisyCellMap((TH2*)ConvertMapRuns2Events(nruns,runNumbers,hBadCellMap[3]), "hMapEvents_ETotalClusterEhigh");
208 // ... or we can also draw the percent for each cell being dead/noisy
209 // DrawDeadNoisyCellPercent(nruns, runNumbers, hBadCellMap[0], "hPercent_NTimesInClusterElow");
210 // DrawDeadNoisyCellPercent(nruns, runNumbers, hBadCellMap[1], "hPercent_NTimesInClusterEhigh");
211 // DrawDeadNoisyCellPercent(nruns, runNumbers, hBadCellMap[2], "hPercent_ETotalClusterElow");
212 // DrawDeadNoisyCellPercent(nruns, runNumbers, hBadCellMap[3], "hPercent_ETotalClusterEhigh");
214 // Our main criteria for analytical finding of bad cells is based on the following.
215 // Note that factor distributions for NTimesInCluster and ETotalCluster are very
216 // similar, both at low and at high energy. Note also that we can say the same
217 // for dead/noisy cell visualizations above: they are very similar. This suggests
218 // that NTimesInCluster and ETotalCluster give the same results. The criteria
219 // NTimesInCluster, however, is more calibration-independent (though energy thresholds
220 // are still calibration-dependent) and thus is more reliable and clear. Thus, we
221 // limit ourselves to NTimesInClusterElow/NTimesInClusterEhigh criteria.
222 // Now, you could probably noted that NTimesInClusterEhigh give more dead
223 // cells than that at low energy. This is an expected statistical effect: we have
224 // much smaller number of clusters at high energy. Consequently, we will not use dead
225 // cell candidates at high energy.
227 // Finally, we combine candidates from low/high energies and produce one TH2
228 // histogram which is the primary source of our analytical results.
230 TH2* hBadCellMapPrimary = CombineDeadNoisyElowEhigh(hBadCellMap[0], hBadCellMap[1]);
232 // Note for PHOS: if you are not happy with NTimesInClusterEhigh results (due to a lack of statistics)
233 // uncomment this line:
234 // hBadCellMapPrimary = hBadCellMap[0];
236 // Draw everything for combined
237 DrawDeadNoisyCellMap(hBadCellMapPrimary, "Primary_hMapRuns");
238 DrawDeadNoisyCellMap((TH2*)ConvertMapRuns2Events(nruns,runNumbers,hBadCellMapPrimary), "Primary_hMapEvents");
239 DrawDeadNoisyCellPercent(nruns, runNumbers, hBadCellMapPrimary, "Primary_hPercent");
241 // Print full information on cells which are dead/noisy;
242 // kTRUE -- print also the percentage % dead/% noisy
243 PrintDeadNoisyCells(hBadCellMapPrimary, 0.9, 1.); // in 90-100% of runs
244 PrintDeadNoisyCells(hBadCellMapPrimary, 0.5, 0.9, kTRUE); // in 50-90% of runs
245 // PrintDeadNoisyCells(hBadCellMapPrimary, 0.3, 0.5); // in 30-50% of runs
247 // visualize dead/noisy cell map for EMCAL/PHOS; requires aliroot
248 // DrawOccupancy(nruns, runNumbers, hBadCellMapPrimary, "hDeadNoisyCellsOccupancy");
250 // EMCAL: print full information on missing/noisy parts (e.g. RCUs); requires aliroot
251 // PrintEMCALProblematicBlocks(nruns, runNumbers, hBadCellMapPrimary);
254 /* RUNS QUALITY SECTION: CLUSTER AVERAGES + PI0 AVERAGES
258 // First, draw cluster averages per run;
259 // 1 = minimum number of cells in a cluster;
260 // 0.3GeV = minimum cluster energy;
261 // -1 = maximum cluster energy = infinity (in fact, 20GeV in the default configuration).
262 DrawClusterAveragesPerRun(nruns, runNumbers, 1, 0.3, -1);
264 // Second, draw the same cluster averages per run, but corrected for detector acceptance
265 DrawClusterAveragesPerRun(nruns, runNumbers, 1, 0.3, -1, hBadCellMapPrimary);
267 // Draw random slices of the pi0 peak in one supermodule and in whole detector
268 DrawPi0Slice(runNumbers[gRandom->Integer(nruns)], 1 + gRandom->Integer(3));
269 DrawPi0Slice(runNumbers[gRandom->Integer(nruns)], -1);
271 // Draw number of pi0s per event, pi0 mass and width
272 DrawPi0Averages(nruns, runNumbers);
274 // Draw pi0 values per event with pi0 photons both in the same supermodule
275 // DrawPi0Averages(nruns, runNumbers, kTRUE);
277 // Draw pi0 values per event with pi0 photons both in the same supermodule
278 // + correct for supermodule acceptance (should not be treated to be 100% reliable!)
279 // DrawPi0Averages(nruns, runNumbers, kTRUE, hBadCellMap[0]);
281 // Draw pi0 distributions (helps to take decision on bad runs);
282 // first argument -- suffix from canvas title;
283 // second argument -- number of bins in distributions
284 // DrawPi0Distributions("", 100);
285 // DrawPi0Distributions("SM1", 100)
286 // DrawPi0Distributions("_sameSM", 100);
287 // DrawPi0Distributions("_sameSM_corr4accep", 100);
294 * Lazy/curious boundary is here! -------------------------
295 * Do not pass it if you do not fill curious enough! ;)
297 * The shape analysis below belongs to the main bad cell searching criteria.
298 * However, it may require some extra work in order to make it usable.
300 * The idea behind is simple: fit the energy spectrum of a cell with
301 * the function A*exp(-B*x)/x^2 (which proved to be a very good fit),
302 * draw parameters A, B and chi2/ndf distributions and notice
303 * cells which deviate too much from the averages.
305 * Huge statistics (>20M events) is necessary for this to work reliably.
310 // The analysis below defines a cell as being bad simply if it is outside
311 // of some good region for any of the parameters A, B, chi2/ndf.
312 // The regions are depicted by vertical orange lines. The problem is that these
313 // regions are usually not automatically determined correctly.
315 // The function syntax is the following:
316 // 0.1-1.0 GeV -- fitting region;
317 // hCellAmplitude -- the histogram which to take for processing,
318 // the other possible choice is hCellAmplitudeNonLocMax;
319 // 1000 -- number of bins is distributions.
321 // The next three groups of parameters are:
322 // <text label> <maximum value on distribution> <left edge of the good region> <right edge of the good region>
324 // It is left/right edges which usually require manual tuning.
325 // -1 means to initialize a parameter automatically.
327 TestCellShapes(0.25, 1., "hCellAmplitude", 1000,
328 // maxval / left margin / right margin
331 -1, -1,-1); // fit chi2/ndf
334 /* DISTRIBUTION ANALYSIS
336 * Test for bad cells by plotting a distribution among cells.
338 * It is especially useful in searching for miscalibrated cells.
339 * The function parameters are similar to parameters from shape analysis section.
342 TestCellsTH1(nruns, runNumbers, "hCellLocMaxNTimesInClusterElow",
343 "Number of times cell was local maximum, low energy, per event", "Entries",
344 400, -1,-1,-1); // nbins / maxval / left margin / right margin
346 TestCellsTH1(nruns, runNumbers, "hCellLocMaxNTimesInClusterEhigh",
347 "Number of times cell was local maximum, high energy, per event", "Entries",
348 400, -1,-1,-1); // nbins / maxval / left margin / right margin
350 // TestCellsTH1(nruns, runNumbers, "hCellLocMaxETotalClusterElow",
351 // "Total cluster energy for local maximum cell, low energy, per event", "Entries",
352 // 400, -1,-1,-1); // nbins / maxval / left margin / right margin
354 // TestCellsTH1(nruns, runNumbers, "hCellLocMaxETotalClusterEhigh",
355 // "Total cluster energy for local maximum cell, high energy, per event", "Entries",
356 // 400, -1,-1,-1); // nbins / maxval / left margin / right margin
358 // Three more tests for bad cells:
359 // 1) total deposited energy;
360 // 2) total number of entries;
361 // 3) average energy = [total deposited energy]/[total number of entries].
363 TestCellEandN(0.1, "hCellAmplitude", 1000,
364 // maxval / left margin / right margin
367 -1,-1,-1); // cell E/N
369 // the same at high energies
370 TestCellEandN(0.1, "hCellAmplitudeEhigh", 1000,
371 // maxval / left margin / right margin
374 -1,-1,-1); // cell E/N
377 // The same as above, but for not a local maximum cells; require more statistics
378 // TestCellsTH1(nruns, runNumbers, "hCellNonLocMaxNTimesInClusterElow",
379 // "Number of times cell wasn't local maximum, low energy, per event", "Entries",
380 // 400, -1,-1,-1); // nbins / maxval / left margin / right margin
382 // TestCellsTH1(nruns, runNumbers, "hCellNonLocMaxNTimesInClusterEhigh",
383 // "Number of times cell wasn't local maximum, high energy, per event", "Entries",
384 // 400, -1,-1,-1); // nbins / maxval / left margin / right margin
386 // TestCellsTH1(nruns, runNumbers, "hCellNonLocMaxETotalClusterElow",
387 // "Total cluster energy for not local maximum cell, low energy, per event", "Entries",
388 // 400, -1,-1,-1); // nbins / maxval / left margin / right margin
390 // TestCellsTH1(nruns, runNumbers, "hCellNonLocMaxETotalClusterEhigh",
391 // "Total cluster energy for not local maximum cell, high energy, per event", "Entries",
392 // 400, -1,-1,-1); // nbins / maxval / left margin / right margin
394 // TestCellEandN(0.1, "hCellAmplitudeNonLocMax", 1000,
395 // // maxval / left margin / right margin
396 // -1,-1,-1, // cell E
397 // -1,-1,-1, // cell N
398 // -1,-1,-1); // cell E/N
400 // TestCellEandN(0.1, "hCellAmplitudeEhighNonLocMax", 1000,
401 // // maxval / left margin / right margin
402 // -1,-1,-1, // cell E
403 // -1,-1,-1, // cell N
404 // -1,-1,-1); // cell E/N
414 /* BAD CELLS SEARCHING FUNCTIONS
418 //_________________________________________________________________________
419 TH2** FindDeadNoisyCellsPerRun(const Int_t nruns, Int_t runNumbers[],
420 Double_t factorDead = 0.01, Double_t factorNoisy = 4.,
421 Int_t fnbins = 200, Double_t fxmax = 10.)
423 // Return bad cell candidate maps for the four criteria;
424 // X axis -- cell absId, Y axis -- run index, content: 0=not bad, 1=noisy, -1=dead.
426 // For each run and each cell calculate factor = [cell value]/[averate over cells],
427 // mark cell as noisy if factor >= factorNoisy, mark cell as dead if factor <= factorDead.
429 // fnbins, fxmax -- number of bins and X axis maximum value for factor distributions.
431 // Four criteria in order:
432 char *hname[4] = {"hCellLocMaxNTimesInClusterElow", "hCellLocMaxNTimesInClusterEhigh",
433 "hCellLocMaxETotalClusterElow", "hCellLocMaxETotalClusterEhigh"};
435 TH1* hFactorDistr[4];
436 TH2** hBadCellMap = new TH2*[4];
438 // take one histogram to get binning parameters
439 TH1* one = (TH1*) gFile->Get(Form("run%i_%s",runNumbers[0],hname[0]));
440 Int_t ncells = one->GetNbinsX();
441 Double_t amin = one->GetXaxis()->GetXmin();
442 Double_t amax = one->GetXaxis()->GetXmax();
444 // find dead/noisy cell candidates
445 for (Int_t k = 0; k < 4; k++) {
446 hBadCellMap[k] = new TH2C(Form("hBadCellMap_%s_fdead=%.3f_fnoisy=%.1f",hname[k],factorDead,factorNoisy),
447 Form("Dead/noisy cell map (%s, fdead=%.3f, fnoisy=%.1f)",hname[k],factorDead,factorNoisy),
448 ncells,amin,amax, nruns,0,nruns);
449 hBadCellMap[k]->SetXTitle("AbsId");
450 hBadCellMap[k]->SetYTitle("Run index");
451 hBadCellMap[k]->SetTitleOffset(1.7,"Y");
453 hFactorDistr[k] = new TH1F(Form("hFactorDistr_%s_fdead=%.3f_fnoisy=%.1f",
454 hname[k],factorDead,factorNoisy), "", fnbins,0,fxmax);
455 hFactorDistr[k]->SetTitle("Factor distributions in all runs");
456 hFactorDistr[k]->SetXTitle("Factor");
457 hFactorDistr[k]->SetYTitle("Entries");
460 for (Int_t ri = 0; ri < nruns; ri++) {
461 TH1* one = (TH1*) gFile->Get(Form("run%i_%s", runNumbers[ri], hname[k]));
464 Double_t av = 0; // average
465 Int_t count = 0; // counted cells
466 for (Int_t c = 1; c <= ncells; c++) {
467 // do not count cells with zero content
468 if (one->GetBinContent(c) == 0) continue;
469 // do not count cells listed in the parameters section in the beginning of the macro
470 if (IsCellMarked4Exclude(one->GetBinLowEdge(c))) continue;
472 av += one->GetBinContent(c);
475 // division by zero checks
477 Warning("FindDeadNoisyCellsPerRun", Form("No cells counted at ri=%i",ri));
483 Warning("FindDeadNoisyCellsPerRun", Form("Average is zero at ri=%i",ri));
487 // find dead/noisy candidates
488 for (Int_t c = 1; c <= ncells; c++) {
489 Double_t fac = one->GetBinContent(c)/av;
490 hFactorDistr[k]->Fill(fac);
492 if (fac <= factorDead)
493 hBadCellMap[k]->SetBinContent(c, ri+1, -1);
494 else if (fac >= factorNoisy)
495 hBadCellMap[k]->SetBinContent(c, ri+1, 1);
503 // draw factor distributions ...
504 TCanvas *c1 = new TCanvas("hFactorDistr", "hFactorDistr", 400,400);
505 gPad->SetLeftMargin(0.12);
506 gPad->SetRightMargin(0.08);
508 hFactorDistr[0]->SetLineColor(kBlue);
509 hFactorDistr[1]->SetLineColor(kRed);
510 hFactorDistr[2]->SetLineColor(kGreen);
511 hFactorDistr[3]->SetLineColor(kOrange);
512 hFactorDistr[0]->Draw();
513 hFactorDistr[1]->Draw("same");
514 hFactorDistr[2]->Draw("same");
515 hFactorDistr[3]->Draw("same");
518 TLegend *leg = new TLegend(0.45,0.65,0.90,0.85);
519 leg->AddEntry(hFactorDistr[0], "NTimesInCluster, low energy","l");
520 leg->AddEntry(hFactorDistr[2], "ETotalCluster, low energy","l");
521 leg->AddEntry(hFactorDistr[1], "NTimesInCluster, high energy","l");
522 leg->AddEntry(hFactorDistr[3], "EtotalCluster, high energy","l");
530 //_________________________________________________________________________
531 void PrintCellsToExcludeFromAverages(TH2** hBadCellMap)
533 // Print cells suggested for exclusion from averages calculation
535 Int_t ncells = hBadCellMap[0]->GetNbinsX();
536 Int_t nruns = hBadCellMap[0]->GetNbinsY();
538 Int_t *suggested = new Int_t[ncells];
539 memset(suggested, 0, ncells*sizeof(Int_t));
541 for (Int_t c = 1; c <= ncells; c++)
542 for (Int_t ri = 1; ri <= nruns; ri++)
543 if (hBadCellMap[0]->GetBinContent(c, ri) != 0 || hBadCellMap[2]->GetBinContent(c, ri) != 0 ||
544 hBadCellMap[1]->GetBinContent(c, ri) > 0 || hBadCellMap[3]->GetBinContent(c, ri) > 0) // NOTE: dead not counted
547 printf("Suggested cells to switch off in averages calculations (copai 'n paste!):\n");
548 printf("Int_t excells[] = {");
550 for (Int_t c = 0; c < ncells; c++)
551 if (suggested[c] >= 0.4*nruns) {// 40% of runs threshold
552 printf("%s%i", n == 0 ? "" : ",", c);
555 printf("};\nInt_t nexc = %i;\n\n",n);
558 //_________________________________________________________________________
559 Bool_t IsCellMarked4Exclude(Int_t absId)
561 // Return true if a cell is in excells[] array
563 static TH1* hExclCells = NULL;
565 // one time initialization
567 hExclCells = new TH1F("hExclCells", "", 20000,0,20000);
568 for (Int_t c = 0; c < nexc; c++)
569 hExclCells->SetBinContent(hExclCells->FindBin(excells[c]), 1);
572 return (hExclCells->GetBinContent(hExclCells->FindBin(absId)) > 0 ? kTRUE : kFALSE);
575 //_________________________________________________________________________
576 void DrawDeadNoisyCellMap(TH2* hmap, char* cname)
578 // Visualize dead/noisy cell map;
579 // cname -- canvas name.
581 TCanvas *c1 = new TCanvas(cname, cname, 0,0,600,600);
582 gPad->SetLeftMargin(0.14);
583 gPad->SetRightMargin(0.06);
585 // draw dummy to initialize axis ranges
586 TH2* hDummy = (TH2*) hmap->Clone(Form("hDummy_%s",cname));
590 for (Int_t c = 1; c <= hmap->GetNbinsX(); c++) //cell number
591 for (Int_t ri = 1; ri <= hmap->GetNbinsY(); ri++) { //run index
592 Double_t stat = hmap->GetBinContent(c, ri); // cell status
593 if (stat == 0) continue;
595 Double_t x = hmap->GetXaxis()->GetBinCenter(c);
596 Double_t y1 = hmap->GetYaxis()->GetBinLowEdge(ri);
597 Double_t y2 = hmap->GetYaxis()->GetBinUpEdge(ri);
599 // draw a line; FIXME: what is a better choice?
600 TLine* line = new TLine(x,y1,x,y2);
601 line->SetLineWidth(1);
602 if (stat > 0) line->SetLineColor(kRed); // noisy cell
603 else line->SetLineColor(kBlue); // dead cell
610 //_________________________________________________________________________
611 void DrawDeadNoisyCellPercent(Int_t nruns, Int_t runNumbers[], TH2* hmap, char* cname)
613 // Show percent of runs/events for each cell being problematic;
614 // cname -- canvas name.
616 // binning parameters
617 Int_t ncells = hmap->GetNbinsX();
618 Double_t amin = hmap->GetXaxis()->GetXmin();
619 Double_t amax = hmap->GetXaxis()->GetXmax();
621 // number of times cell was dead/noisy;
622 Int_t *ndeadRuns = new Int_t[ncells];
623 Int_t *nnoisyRuns = new Int_t[ncells];
624 Double_t *ndeadEvents = new Double_t[ncells];
625 Double_t *nnoisyEvents = new Double_t[ncells];
628 for (Int_t c = 0; c < ncells; c++) {
634 for (Int_t ri = 0; ri < nruns; ri++) {
635 Double_t stat = hmap->GetBinContent(c+1, ri+1); // cell status
636 Int_t nevents = GetNumberOfEvents(runNumbers[ri]);
640 nnoisyEvents[c] += nevents;
644 ndeadEvents[c] += nevents;
649 // total number of events
650 Double_t ntotal = GetTotalNumberOfEvents(nruns, runNumbers);
652 TCanvas *c1 = new TCanvas(cname, cname, 0,0,600,600);
653 gPad->SetLeftMargin(0.14);
654 gPad->SetRightMargin(0.06);
656 // draw dummy histogram to initialize canvas
657 TH1* hDummy = new TH1F(Form("hDummy_%s",cname), hmap->GetTitle(), ncells,amin,amax);
658 hDummy->SetAxisRange(0,100, "Y");
659 hDummy->SetXTitle("AbsId");
660 hDummy->SetYTitle("Percent");
669 // draw results, FIXME: is where a better way?
670 for (Int_t c = 0; c < ncells; c++) {
671 Double_t x = hmap->GetXaxis()->GetBinCenter(c+1);
672 Double_t y1 = 100.*ndeadEvents[c]/ntotal;
673 Double_t y2 = 100.*ndeadRuns[c]/nruns;
676 if (ndeadEvents[c] > 0) {
677 line1 = new TLine(x, 0, x, y1);
678 line1->SetLineWidth(1);
679 line1->SetLineColor(kBlue);
684 if (nnoisyEvents[c] > 0) {
685 line2 = new TLine(x, y1, x, y1 + 100.*nnoisyEvents[c]/ntotal);
686 line2->SetLineWidth(1);
687 line2->SetLineColor(kRed);
692 if (ndeadRuns[c] > 0) {
693 line3 = new TLine(x, 0, x, y2);
694 line3->SetLineWidth(1);
695 line3->SetLineStyle(7);
696 line3->SetLineColor(kBlack);
701 if (nnoisyRuns[c] > 0) {
702 line4 = new TLine(x, y2, x, y2 + 100.*nnoisyRuns[c]/nruns);
703 line4->SetLineWidth(1);
704 line4->SetLineStyle(7);
705 line4->SetLineColor(kOrange);
711 TLegend *leg = new TLegend(0.65,0.65,0.92,0.75);
712 if (line1) leg->AddEntry(line1, "% of events, dead","l");
713 if (line2) leg->AddEntry(line2, "% of events, noisy","l");
714 if (line3) leg->AddEntry(line3, "% of runs, dead","l");
715 if (line4) leg->AddEntry(line4, "% of runs, noisy","l");
721 //_________________________________________________________________________
722 TH1* ConvertMapRuns2Events(const Int_t nruns, Int_t runNumbers[], TH1* inhisto)
724 // Returns a histogram in which run index axis is converted to number of events
725 // by making a variable axis bin width.
726 // If inhisto inherits from TH2, convert Y axis; convert X axis otherwise.
729 Double_t *nevents = new Double_t[nruns+1];
731 for (Int_t ri = 0; ri < nruns; ri++)
732 nevents[1+ri] = nevents[ri] + GetNumberOfEvents(runNumbers[ri]);
734 TH1* outh = (TH1*) inhisto->Clone(Form("%s_Events",inhisto->GetName()));
736 if (inhisto->InheritsFrom("TH2")) {
737 outh->GetYaxis()->Set(nruns, nevents);
738 outh->SetYTitle("Number of events");
739 outh->SetTitleOffset(1.7,"Y");
742 outh->GetXaxis()->Set(nruns, nevents);
743 outh->SetXTitle("Number of events");
749 //_________________________________________________________________________
750 TH2* CombineDeadNoisyElowEhigh(TH2* hmapElow, TH2* hmapEhigh)
752 // Combine two maps at low/high energy into one,
753 // do not count dead map at high energy.
754 // NOTE: if cell is dead at Elow and noisy and Ehigh, set content = 2.
756 TH2* hmap_combined = (TH2*) hmapElow->Clone(Form("%s+%s",hmapElow->GetName(),hmapEhigh->GetName()));
758 for (Int_t c = 1; c <= hmapElow->GetNbinsX(); c++)
759 for (Int_t ri = 1; ri <= hmapElow->GetNbinsY(); ri++)
760 if (hmapEhigh->GetBinContent(c, ri) > 0) {
761 hmap_combined->SetBinContent(c, ri, 1);
763 if (hmapElow->GetBinContent(c, ri) < 0)
764 hmap_combined->SetBinContent(c, ri, 2);
767 hmap_combined->SetTitle("Dead/noisy cell map, combined");
769 return hmap_combined;
772 //_________________________________________________________________________
773 void PrintDeadNoisyCells(TH2* hmap, Double_t percent1 = 0.95, Double_t percent2 = 1., Bool_t kPrintPercentage = kFALSE)
775 // Print full information on dead/noisy cells which are present in
776 // percent1-percent2 portion of runs (percent1 excluded, percent2 included).
778 Int_t ncells = hmap->GetNbinsX();
779 Int_t nruns = hmap->GetNbinsY();
783 printf("Dead/noisy cells in >%.1f%% and <=%.1f%% of runs:", 100*percent1, 100*percent2);
785 for (Int_t c = 1; c <= ncells; c++) {
786 Int_t nrdead = 0; // count number of runs for the current cell
789 for (Int_t ri = 1; ri <= nruns; ri++)
790 if (hmap->GetBinContent(c,ri) < 0) nrdead++;
791 else if (hmap->GetBinContent(c,ri) > 0) nrnoisy++;
793 if (nrdead+nrnoisy > percent1*nruns && nrdead+nrnoisy <= percent2*nruns) {
794 printf(" %.0f", hmap->GetBinLowEdge(c));
795 if (kPrintPercentage)
796 printf("(%-2.0f/%-2.0f)", 100.*nrdead/nruns, 100.*nrnoisy/nruns);
801 printf(" (%i total)\n\n", nbad);
804 //_________________________________________________________________________
805 void DrawOccupancy(Int_t nruns, Int_t runNumbers[], TH2* hmap, char* cname)
807 // Draw bad cell map for EMCAL or PHOS;
808 // cname -- canvas name.
811 if (hmap->GetNbinsX() % 1152 == 0)
812 DrawEMCALOccupancy(nruns, runNumbers, hmap, cname);
814 DrawPHOSOccupancy(nruns, runNumbers, hmap, cname);
817 //_________________________________________________________________________
818 void DrawEMCALOccupancy(Int_t nruns, Int_t runNumbers[], TH2* hmap, char* cname)
820 // Draw bad cell map for EMCAL;
821 // cname -- canvas name.
823 Int_t nmods = hmap->GetNbinsX()/1152; // number of supermodules
824 Int_t vsize = ceil(nmods/2.); // vertical size in canvas
825 Int_t lastSector = (nmods-1)/2; // for pad number calculation
827 TCanvas *c1 = new TCanvas(cname, cname, 800,200*vsize);
828 c1->Divide(2, vsize);
830 for (Int_t sm = 0; sm < nmods; sm++)
832 // top left is SM0, bottom right is SM9
835 c1->cd((lastSector - sector)*2 + side + 1);
837 TH2* hSM = new TH2C(Form("hEMCALSM%i_%s",sm,cname), Form("SM%i",sm), 48,0,48, 24,0,24);
838 hSM->SetXTitle("iEta");
839 hSM->SetYTitle("iPhi");
841 // loop over supermodule cells
842 for (Int_t c = 0; c < 1152; c++) {
843 Int_t absId = 1152 * sm + c;
845 for (Int_t ri = 0; ri < nruns; ri++) {
846 if (hmap->GetBinContent(hmap->FindBin(absId,ri)) == 0) continue;
848 Int_t nModule, eta, phi;
849 AbsId2EtaPhi(absId, nModule, eta, phi);
855 } // supermodule loop
860 //_________________________________________________________________________
861 void DrawPHOSOccupancy(Int_t nruns, Int_t runNumbers[], TH2* hmap, char* cname)
863 // Draw bad cell map for PHOS;
864 // cname -- canvas name.
866 Int_t nmods = hmap->GetNbinsX()/3584; // number of supermodules
867 Int_t vsize = nmods; // vertical size in canvas
869 TCanvas *c1 = new TCanvas(cname, cname, 64*5,56*5*vsize);
870 c1->Divide(1, vsize);
872 for (Int_t sm = 1; sm <= nmods; sm++)
876 TH2* hSM = new TH2C(Form("hPHOSSM%i_%s",sm,cname), Form("SM%i",sm), 64,0,64, 56,0,56);
877 hSM->SetXTitle("iEta");
878 hSM->SetYTitle("iPhi");
880 // loop over supermodule cells
881 for (Int_t c = 1; c <= 3584; c++) {
882 Int_t absId = 3584*(sm-1) + c;
884 for (Int_t ri = 0; ri < nruns; ri++) {
885 if (hmap->GetBinContent(hmap->FindBin(absId,ri)) == 0) continue;
887 Int_t nModule, eta, phi;
888 AbsId2EtaPhi(absId, nModule, eta, phi, 1);
894 } // supermodule loop
899 //_________________________________________________________________________
900 void PrintEMCALProblematicBlocks(Int_t nruns, Int_t runNumbers[], TH2* hmap)
902 // Print, on a per run basis, complete information about EMCAL missing
903 // (or dead or noisy) blocks. Missing/noisy EMCAL atomic part is a 2x2
904 // block (288 blocks per supermodule).
906 // number of supermodules
907 Int_t nmods = hmap->GetNbinsX()/1152;
909 Printf("Problematic (missing/dead/noisy) 2x2 block numbers in EMCAL:");
912 for (Int_t ri = 0; ri < nruns; ri++) {
913 Printf("Run %i:", runNumbers[ri]);
916 for (Int_t sm = 0; sm < nmods; sm++) {
918 // will be filled with the number of missing cells (from 0 to 4)
920 memset(blk2x2, 0, 288*sizeof(Int_t));
922 for (Int_t c = 0; c < 1152; c++) {
923 Int_t absId = 1152 * sm + c;
925 // select problematic cells only
926 if (hmap->GetBinContent(hmap->FindBin(absId,ri)) == 0) continue;
928 Int_t nModule, eta, phi;
929 AbsId2EtaPhi(c, nModule, eta, phi);
933 // calculate number of bad blocks
935 for (Int_t b = 0; b < 288; b++)
936 if (blk2x2[b] == 4) nbad2x2++;
939 if (nbad2x2 == 0) continue;
941 printf(" SM%i:", sm);
944 if (nbad2x2 == 288) {
945 printf(" missing the whole supermodule!\n");
950 if (nbad2x2 >= 144) {
955 for (Int_t b = 0; b < 288; b++)
956 if (blk2x2[b] == 4) nRCU[GetEMCALRCUNumber(b)]++;
958 if (nRCU[0] == 144) {
960 for (Int_t b = 0; b < 288; b++)
961 if (blk2x2[b] == 4 && GetEMCALRCUNumber(b) == 0) blk2x2[b] = 0;
964 if (nRCU[1] == 144) {
966 for (Int_t b = 0; b < 288; b++)
967 if (blk2x2[b] == 4 && GetEMCALRCUNumber(b) == 1) blk2x2[b] = 0;
975 for (Int_t b = 0; b < 288; b++)
976 if (blk2x2[b] == 4) printf(" %i", b);
977 printf(" (%i)", nbad2x2);
982 } // supermodule loop
988 //_________________________________________________________________________
989 Int_t GetEMCALRCUNumber(Int_t nModule)
991 // Returns RCU number for a 2x2 block in EMCAL;
992 // nModule -- block number (0-287).
994 static TH1* hRCUs = NULL;
996 // one-time initialization
998 hRCUs = new TH1F("hRCU1", "", 288,0,288);
1000 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,
1001 92,93,94,95,104,105,106,107,116,117,118,119,128,129,130,131,140,141,142,143,148,149,150,151,
1002 152,153,154,155,160,161,162,163,164,165,166,167,172,173,174,175,176,177,178,179,184,185,186,
1003 187,188,189,190,191,196,197,198,199,200,201,202,203,208,209,210,211,212,213,214,215,220,221,
1004 222,223,224,225,226,227,232,233,234,235,236,237,238,239,244,245,246,247,248,249,250,251,256,
1005 257,258,259,260,261,262,263,268,269,270,271,272,273,274,275,280,281,282,283,284,285,286,287};
1007 for (Int_t i = 0; i < 144; i++) hRCUs->SetBinContent(RCU1[i]+1, 1);
1010 return hRCUs->GetBinContent(nModule+1);
1013 //_________________________________________________________________________
1014 void AbsId2EtaPhi(Int_t absId, Int_t &nModule, Int_t &eta, Int_t &phi, Int_t det = 0)
1016 // Converts cell absId --> (sm,eta,phi);
1018 // nModule -- 2x2 block number for EMCAL;
1019 // det -- detector: 0/EMCAL, 1/PHOS.
1023 AliEMCALGeometry *geomEMCAL = AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
1025 Int_t nSupMod, nIphi, nIeta;
1026 geomEMCAL->GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
1027 geomEMCAL->GetCellPhiEtaIndexInSModule(nSupMod, nModule, nIphi, nIeta, phi, eta);
1032 else if (det == 1) {
1033 AliPHOSGeometry *geomPHOS = AliPHOSGeometry::GetInstance("IHEP");
1036 geomPHOS->AbsToRelNumbering(absId, relid);
1047 Error("AbsId2EtaPhi", "Wrong detector");
1051 //_________________________________________________________________________
1052 void TestCellsTH1(Int_t nruns, Int_t runNumbers[], char *hname,
1053 char* htitle = "", char* hytitle = "",
1054 Int_t dnbins = 200, Double_t dmaxval = -1, Double_t goodmin = -1, Double_t goodmax = -1)
1056 // Test for bad cells by plotting a distribution of a TH1 histogram.
1057 // The histogram is obtained as a sum over runs of TH1 per run.
1059 // hname -- histogram name to process;
1060 // htitle, hytitle -- histogram title and Y axis title;
1061 // dnbins -- number of bins in distribution;
1062 // dmaxval -- X axis maximum in distribution.
1063 // goodmin,goodmax -- the region outside which a cell is considered as bad.
1065 // -1 value for maxval/goodmin/goodmax -- process a variable automatically.
1067 // initialize histogram
1068 TH1* histo = (TH1*) gFile->Get(Form("run%i_%s", runNumbers[0], hname));
1069 histo->SetName(hname);
1070 histo->SetTitle(htitle);
1071 histo->SetXTitle("AbsId");
1072 histo->SetYTitle(hytitle);
1075 for (Int_t i = 1; i < nruns; i++) {
1076 TH1* h = (TH1*) gFile->Get(Form("run%i_%s", runNumbers[i], hname));
1081 histo->Scale(1./(Double_t)GetTotalNumberOfEvents(nruns, runNumbers));
1082 Process1(histo, Form("TestCellsTH1_%s",hname), dnbins, dmaxval, goodmin, goodmax);
1085 //_________________________________________________________________________
1086 void TestCellEandN(Double_t Emin = 0.1, char* hname = "hCellAmplitude", Int_t dnbins = 200,
1087 Double_t maxval1 = -1, Double_t goodmin1 = -1, Double_t goodmax1 = -1,
1088 Double_t maxval2 = -1, Double_t goodmin2 = -1, Double_t goodmax2 = -1,
1089 Double_t maxval3 = -1, Double_t goodmin3 = -1, Double_t goodmax3 = -1)
1091 // Three more tests for bad cells:
1092 // 1) total deposited energy;
1093 // 2) total number of entries;
1094 // 3) average energy = [total deposited energy]/[total number of entries].
1096 // Based on summary histograms. Possible choises:
1097 // hCellAmplitude, hCellAmplitudeEhigh, hCellAmplitudeNonLocMax, hCellAmplitudeEhighNonLocMax
1099 // Emin -- minimum cell amplitude to count;
1100 // hname -- name (in file) of TH2 histogram with cell amplitudes;
1101 // dnbins -- number of bins in distributions;
1102 // maxval[123] -- maximum values on distributions for the criteria 1),2),3) respectively;
1103 // goodmin[123],goodmax[123] -- the regions on distributions outside those a cell is considered as bad.
1105 // input; X axis -- absId numbers
1106 TH2 *hCellAmplitude = (TH2*) gFile->Get(hname);
1108 // binning parameters
1109 Int_t ncells = hCellAmplitude->GetNbinsX();
1110 Double_t amin = hCellAmplitude->GetXaxis()->GetXmin();
1111 Double_t amax = hCellAmplitude->GetXaxis()->GetXmax();
1113 TH1* hCellEtotal = new TH1F(Form("%s_hCellEtotal_E%.2f",hname,Emin),
1114 Form("Total deposited energy, E > %.2f GeV",Emin), ncells,amin,amax);
1115 hCellEtotal->SetXTitle("AbsId");
1116 hCellEtotal->SetYTitle("Energy, GeV");
1118 TH1F *hCellNtotal = new TH1F(Form("%s_hCellNtotal_E%.2f",hname,Emin),
1119 Form("Total number of entries, E > %.2f GeV",Emin), ncells,amin,amax);
1120 hCellNtotal->SetXTitle("AbsId");
1121 hCellNtotal->SetYTitle("Entries");
1123 TH1F *hCellEtoNtotal = new TH1F(Form("%s_hCellEtoNtotal_E%.2f",hname,Emin),
1124 Form("Average energy per hit, E > %.2f GeV",Emin), ncells,amin,amax);
1125 hCellEtoNtotal->SetXTitle("AbsId");
1126 hCellEtoNtotal->SetYTitle("Energy, GeV");
1129 for (Int_t c = 1; c <= ncells; c++) {
1133 for (Int_t i = 1; i <= hCellAmplitude->GetNbinsY(); i++) {
1134 Double_t E = hCellAmplitude->GetYaxis()->GetBinCenter(i);
1135 Double_t N = hCellAmplitude->GetBinContent(c, i);
1136 if (E < Emin) continue;
1141 hCellEtotal->SetBinContent(c, Esum);
1142 hCellNtotal->SetBinContent(c, Nsum);
1144 if (Nsum > 0.5) // number of entries >= 1
1145 hCellEtoNtotal->SetBinContent(c, Esum/Nsum);
1148 delete hCellAmplitude;
1150 Process1(hCellEtotal, Form("%s_CellE", hname), dnbins, maxval1, goodmin1, goodmax1);
1151 Process1(hCellNtotal, Form("%s_CellN", hname), dnbins, maxval2, goodmin2, goodmax2);
1152 Process1(hCellEtoNtotal, Form("%s_CellE/N", hname), dnbins, maxval3, goodmin3, goodmax3);
1155 //_________________________________________________________________________
1156 void TestCellShapes(Double_t fitEmin, Double_t fitEmax, char* hname = "hCellAmplitude", Int_t dnbins = 1000,
1157 Double_t maxval1 = -1, Double_t goodmin1 = -1, Double_t goodmax1 = -1,
1158 Double_t maxval2 = -1, Double_t goodmin2 = -1, Double_t goodmax2 = -1,
1159 Double_t maxval3 = -1, Double_t goodmin3 = -1, Double_t goodmax3 = -1)
1161 // Test cells shape using fit function f(x)=A*exp(-B*x)/x^2.
1162 // Produce values per cell + distributions for A, B and chi2/ndf parameters.
1164 // fitEmin, fitEmax -- fit range;
1165 // hname -- name (in file) of TH2 histogram with cell amplitudes;
1166 // dnbins -- number of bins in distributions;
1167 // maxval[123] -- maximum values on distributions for the criteria A, B and chi2/ndf respectively;
1168 // goodmin[123],goodmax[123] -- the regions on distributions outside those a cell is considered as bad.
1170 // -1 value for maxval/goodmin/goodmax -- process a variable automatically.
1172 // Note: numbers are optimized for EMCAL.
1173 // TODO: check for PHOS
1175 // input; X axis -- absId numbers
1176 TH2 *hCellAmplitude = (TH2*) gFile->Get(hname);
1178 // binning parameters
1179 Int_t ncells = hCellAmplitude->GetNbinsX();
1180 Double_t amin = hCellAmplitude->GetXaxis()->GetXmin();
1181 Double_t amax = hCellAmplitude->GetXaxis()->GetXmax();
1183 // initialize histograms
1184 TH1 *hFitA = new TH1F(Form("hFitA_%s",hname),"Fit A value", ncells,amin,amax);
1185 hFitA->SetXTitle("AbsId");
1186 hFitA->SetYTitle("A");
1188 TH1 *hFitB = new TH1F(Form("hFitB_%s",hname),"Fit B value", ncells,amin,amax);
1189 hFitB->SetXTitle("AbsId");
1190 hFitB->SetYTitle("B");
1192 TH1 *hFitChi2Ndf = new TH1F(Form("hFitChi2Ndf_%s",hname),"Fit #chi^{2}/ndf value", ncells,amin,amax);
1193 hFitChi2Ndf->SetXTitle("AbsId");
1194 hFitChi2Ndf->SetYTitle("#chi^{2}/ndf");
1196 // total number of events; to estimate A value
1197 TH1* hNEventsProcessedPerRun = (TH1*) gFile->Get("hNEventsProcessedPerRun");
1198 Double_t ntotal = hNEventsProcessedPerRun->Integral(1, hNEventsProcessedPerRun->GetNbinsX());
1201 TF1 *fit = new TF1("fit", "[0]*exp(-[1]*x)/x^2");
1202 fit->SetParLimits(0, ntotal*1e-8,ntotal*1e-4);
1203 fit->SetParLimits(1, 0.,30.);
1204 fit->SetParameter(0, ntotal*1e-6);
1205 fit->SetParameter(1, 1.5);
1207 for (Int_t i = 1; i <= ncells; i++) {
1208 TH1 *hCell = hCellAmplitude->ProjectionY("",i,i);
1209 if (hCell->GetEntries() == 0) continue;
1211 hCell->Fit(fit, "0LQEM", "", fitEmin, fitEmax);
1214 hFitA->SetBinContent(i, fit->GetParameter(0));
1215 hFitB->SetBinContent(i, fit->GetParameter(1));
1216 if (fit->GetNDF() != 0)
1217 hFitChi2Ndf->SetBinContent(i, fit->GetChisquare()/fit->GetNDF());
1220 delete hCellAmplitude;
1222 // automatic parameters, if requested
1223 if (maxval1 < 0) maxval1 = 4e-6 * ntotal;
1224 if (maxval2 < 0) maxval2 = 10.;
1225 if (maxval3 < 0) maxval3 = 15.;
1227 Process1(hFitA, Form("%s_FitA", hname), dnbins, maxval1, goodmin1, goodmax1);
1228 Process1(hFitB, Form("%s_FitB", hname), dnbins, maxval2, goodmin2, goodmax2);
1229 Process1(hFitChi2Ndf, Form("%s_FitChi2ndf", hname), dnbins, maxval3, goodmin3, goodmax3);
1232 //_________________________________________________________________________
1233 void Process1(TH1* inhisto, char* label = "", Int_t dnbins = 200,
1234 Double_t dmaxval = -1, Double_t goodmin = -1, Double_t goodmax = -1)
1236 // Actual distribution analysis for a TH1 histogram:
1237 // 1) create a distribution for the input histogram;
1239 // 3) take a decision about bad cells.
1241 // inhisto -- input histogram;
1242 // label -- text label for stdout;
1243 // dnbins -- number of bins in distribution;
1244 // goodmin,goodmax -- cells outside this region are considered as bad;
1245 // dmaxval -- maximum value on distribution histogram.
1246 // The later is required in cases where a bad cell kills the whole distribution:
1247 // limiting distribution maximum value solves the problem.
1250 dmaxval = inhisto->GetMaximum()*1.01; // 1.01 - to see the last bin
1252 TH1 *distrib = new TH1F(Form("%sDistr",inhisto->GetName()), "", dnbins, inhisto->GetMinimum(), dmaxval);
1253 distrib->SetXTitle(inhisto->GetYaxis()->GetTitle());
1254 distrib->SetYTitle("Entries");
1256 // fill distribution
1257 for (Int_t c = 1; c <= inhisto->GetNbinsX(); c++)
1258 distrib->Fill(inhisto->GetBinContent(c));
1260 // draw histogram + distribution
1261 TCanvas *c1 = new TCanvas(inhisto->GetName(),inhisto->GetName(), 800,400);
1265 gPad->SetLeftMargin(0.14);
1266 gPad->SetRightMargin(0.06);
1268 inhisto->SetTitleOffset(1.7,"Y");
1272 gPad->SetLeftMargin(0.14);
1273 gPad->SetRightMargin(0.10);
1277 // simple way to estimate the left margin for good cells region:
1278 // go from left to right and find the first bin with content 2,
1279 // then go from this bin right to left while bin content is nonzero
1281 goodmin = distrib->GetXaxis()->GetXmin();
1283 for (Int_t i = 1; i <= distrib->GetNbinsX(); i++)
1284 if (distrib->GetBinContent(i) == 2) {
1285 while (i > 1 && distrib->GetBinContent(i-1) > 0) i--;
1286 goodmin = distrib->GetBinLowEdge(i);
1291 // the same automatic algorithm as above, but reflected
1293 goodmax = distrib->GetXaxis()->GetXmax();
1295 for (Int_t i = distrib->GetNbinsX(); i >= 1; i--)
1296 if (distrib->GetBinContent(i) == 2) {
1297 while (i < distrib->GetNbinsX() && distrib->GetBinContent(i+1) > 0) i++;
1298 goodmax = distrib->GetXaxis()->GetBinUpEdge(i);
1304 TLine *lline = new TLine(goodmin, 0, goodmin, distrib->GetMaximum());
1305 lline->SetLineColor(kOrange);
1306 lline->SetLineStyle(7);
1309 TLine *rline = new TLine(goodmax, 0, goodmax, distrib->GetMaximum());
1310 rline->SetLineColor(kOrange);
1311 rline->SetLineStyle(7);
1315 TLegend *leg = new TLegend(0.60,0.82,0.9,0.88);
1316 leg->AddEntry(lline, "Good region boundary","l");
1322 printf("Bad cells by criterum \"%s\":", label);
1324 // print bad cell numbers (outside goodmin,goodmax region)
1326 for (Int_t c = 1; c <= inhisto->GetNbinsX(); c++)
1327 if (inhisto->GetBinContent(c) < goodmin || inhisto->GetBinContent(c) > goodmax) {
1328 printf(" %.0f", inhisto->GetBinLowEdge(c));
1332 printf(" (%i total)\n\n", ntot);
1335 //_________________________________________________________________________
1336 void DrawCell(Int_t absId, Double_t fitEmin = 0.3, Double_t fitEmax = 1.,
1337 Double_t Emax = -1, char* hname = "hCellAmplitude")
1339 // Draw one cell spectrum with a fit.
1341 // fitEmin, fitEmax -- fit range;
1342 // Emax -- maximum value on X axis to show (-1 = no limit);
1343 // hname -- TH2 histogram name to process, possible values:
1344 // "hCellAmplitude", "hCellAmplitudeEHigh", "hCellAmplitudeNonLocMax", "hCellAmplitudeEhighNonLocMax".
1346 // input; X axis -- absId numbers
1347 TH2* hCellAmplitude = (TH2*) gFile->Get(hname);
1349 Int_t bin = hCellAmplitude->GetXaxis()->FindBin(absId);
1350 TH1* hCell = hCellAmplitude->ProjectionY(Form("hCell%i_%s",absId,hname),bin,bin);
1351 hCell->SetXTitle("Energy, GeV");
1352 hCell->SetYTitle("Entries");
1353 hCell->SetTitle(Form("Cell %i, %s", absId, hname));
1354 if (Emax > 0) hCell->SetAxisRange(0, Emax);
1357 TCanvas *c1 = new TCanvas(Form("hCell%i_%s",absId,hname), Form("hCell%i_%s",absId,hname), 400,400);
1358 gPad->SetLeftMargin(0.12);
1359 gPad->SetRightMargin(0.08);
1360 gPad->SetBottomMargin(0.12);
1364 // total number of events
1365 TH1* hNEventsProcessedPerRun = (TH1*) gFile->Get("hNEventsProcessedPerRun");
1366 Double_t ntotal = hNEventsProcessedPerRun->Integral(1, hNEventsProcessedPerRun->GetNbinsX());
1369 TF1 *fit = new TF1("fit", "[0]*exp(-[1]*x)/x^2");
1370 fit->SetLineColor(kRed);
1371 fit->SetLineWidth(2);
1372 fit->SetParName(0, "A");
1373 fit->SetParName(1, "B");
1374 fit->SetParLimits(0, ntotal*1e-8,ntotal*1e-4);
1375 fit->SetParLimits(1, 0.,30.);
1376 fit->SetParameter(0, ntotal*1e-6);
1377 fit->SetParameter(1, 1.);
1378 hCell->Fit(fit,"LQEM", "", fitEmin, fitEmax);
1381 TLegend *leg = new TLegend(0.5,0.75,0.9,0.8);
1382 leg->AddEntry(fit, "A*exp(-B*x)/x^{2}","l");
1388 //_________________________________________________________________________
1389 void DrawCellTime(Int_t absId)
1391 // Draw one cell time spectrum
1393 // input; X axis -- absId numbers
1394 TH2* hCellTime = (TH2*) gFile->Get("hCellTime");
1396 Int_t bin = hCellTime->GetXaxis()->FindBin(absId);
1397 TH1* hCell = hCellTime->ProjectionY(Form("hCellTime%i",absId),bin,bin);
1398 hCell->SetXTitle("Time, s");
1399 hCell->SetYTitle("Entries");
1400 hCell->SetTitle(Form("Cell %i time", absId));
1403 TCanvas *c1 = new TCanvas(Form("hCellTime%i",absId), Form("hCellTime%i",absId), 400,400);
1404 gPad->SetLeftMargin(0.12);
1405 gPad->SetRightMargin(0.10);
1406 gPad->SetBottomMargin(0.12);
1414 /* RUN AVERAGES AND RELATED FUNCTIONS
1418 //_________________________________________________________________________
1419 void DrawClusterAveragesPerRun(Int_t nruns, Int_t runNumbers[], Int_t ncellsMin = 1,
1420 Double_t eclusMin = 0.3, Double_t eclusMax = -1,
1423 // Draws cluster averages per run vs run index and number of events.
1424 // NOTE: the averages are "smoothed" a little due to a finite bin width.
1426 // ncellsMin -- minimum number of cells in cluster (>= 1);
1427 // eclusMin,eclusMax -- minimum/maximum cluster energy cut
1428 // (eclusMax = -1 means infinity);
1429 // hmap -- dead/noisy cell map, which is used for acceptance calculation due
1430 // to switched off detector parts. Acceptance is used to correct the average
1431 // number of clusters per event.
1434 TString s(Form("_NC%i_Emin=%.2fGeV",ncellsMin,eclusMin));
1435 if (eclusMax > 0) s += TString(Form("_Emax=%.2fGeV",eclusMax)).Data();
1436 if (hmap) s += "_corr4accept";
1437 char *suff = s.Data();
1439 if (eclusMax < 0) eclusMax = 1e+5;
1441 // supermodule region
1444 while (SM1 <= SM2 && !gFile->Get(Form("run%i_hNCellsInClusterSM%i",runNumbers[0],SM1)) ) SM1++;
1445 while (SM2 >= SM1 && !gFile->Get(Form("run%i_hNCellsInClusterSM%i",runNumbers[0],SM2)) ) SM2--;
1447 // initialize histograms
1448 hAvECluster = new TH1F(Form("hAvECluster%s",suff), "Average cluster energy", nruns,0,nruns);
1449 hAvECluster->SetXTitle("Run index");
1450 hAvECluster->SetYTitle("Energy, GeV");
1452 hAvNCluster = new TH1F(Form("hAvNCluster%s",suff), "Average number of clusters per event", nruns,0,nruns);
1453 hAvNCluster->SetXTitle("Run index");
1454 hAvNCluster->SetYTitle("Number of clusters");
1456 hAvNCellsInCluster = new TH1F(Form("hAvNCellsInCluster%s",suff), "Average number of cells in cluster", nruns,0,nruns);
1457 hAvNCellsInCluster->SetXTitle("Run index");
1458 hAvNCellsInCluster->SetYTitle("Number of cells");
1460 // initialize per SM histograms
1461 TH1* hAvEClusterSM[10];
1462 TH1* hAvNClusterSM[10];
1463 TH1* hAvNCellsInClusterSM[10];
1465 for (Int_t sm = SM1; sm <= SM2; sm++) {
1466 hAvEClusterSM[sm] = new TH1F(Form("hAvEClusterSM%i%s",sm,suff),"", nruns,0,nruns);
1467 hAvNClusterSM[sm] = new TH1F(Form("hAvNClusterSM%i%s",sm,suff),"", nruns,0,nruns);
1468 hAvNCellsInClusterSM[sm] = new TH1F(Form("hAvNCellsInClusterSM%i%s",sm,suff),"", nruns,0,nruns);
1471 // fill all the histograms per run index
1472 for (Int_t ri = 0; ri < nruns; ri++)
1474 Int_t nevents = GetNumberOfEvents(runNumbers[ri]);
1476 // number of switched off supermodules
1479 Double_t Eclus_total = 0; // total cluster energy
1480 Double_t Nclus_total = 0; // total number of clusters
1481 Double_t Ncells_total = 0; // total number of cells
1484 for (Int_t sm = SM1; sm <= SM2; sm++) {
1485 TH2* hNCellsInClusterSM = (TH2*) gFile->Get(Form("run%i_hNCellsInClusterSM%i",runNumbers[ri],sm));
1487 // the same as above, but per supermodule
1488 Double_t Eclus_totalSM = 0;
1489 Double_t Nclus_totalSM = 0;
1490 Double_t Ncells_totalSM = 0;
1492 // X axis -- cluster energy, Y axis -- number of cells
1493 for (Int_t x = 1; x <= hNCellsInClusterSM->GetNbinsX(); x++)
1494 for (Int_t y = 1+ncellsMin; y <= hNCellsInClusterSM->GetNbinsY(); y++) {//NOTE: bin 1 correspond to ncellsMin=0
1495 Double_t Eclus = hNCellsInClusterSM->GetXaxis()->GetBinCenter(x);
1496 Double_t Ncells = hNCellsInClusterSM->GetYaxis()->GetBinLowEdge(y);
1497 Double_t Nclus = hNCellsInClusterSM->GetBinContent(x,y);
1499 // cut on cluster energy
1500 if (Eclus < eclusMin || Eclus > eclusMax) continue;
1502 Eclus_totalSM += Eclus * Nclus;
1503 Nclus_totalSM += Nclus;
1504 Ncells_totalSM += Ncells * Nclus;
1507 delete hNCellsInClusterSM;
1509 // correct for acceptance
1511 Double_t accep = GetAcceptance(sm, hmap, ri);
1513 Eclus_totalSM /= accep;
1514 Nclus_totalSM /= accep;
1515 Ncells_totalSM /= accep;
1520 Eclus_total += Eclus_totalSM;
1521 Nclus_total += Nclus_totalSM;
1522 Ncells_total += Ncells_totalSM;
1524 hAvNClusterSM[sm]->SetBinContent(ri+1, Nclus_totalSM/nevents);
1525 if (Nclus_totalSM > 0) hAvEClusterSM[sm]->SetBinContent(ri+1, Eclus_totalSM/Nclus_totalSM);
1526 if (Nclus_totalSM > 0) hAvNCellsInClusterSM[sm]->SetBinContent(ri+1, Ncells_totalSM/Nclus_totalSM);
1527 } // supermodule loop
1529 hAvNCluster->SetBinContent(ri+1, Nclus_total/nevents/(SM2-SM1+1-noSM));
1530 if (Nclus_total > 0) hAvECluster->SetBinContent(ri+1, Eclus_total/Nclus_total);
1531 if (Nclus_total > 0) hAvNCellsInCluster->SetBinContent(ri+1, Ncells_total/Nclus_total);
1535 /* Draw results vs run index
1538 TCanvas *c1 = new TCanvas(Form("ClusterAveragesRuns%s",suff),
1539 Form("ClusterAveragesRuns%s",suff), 999,333);
1542 // average cluster energy
1544 gPad->SetLeftMargin(0.14);
1545 gPad->SetRightMargin(0.06);
1546 TLegend *leg = new TLegend(0.65,0.15,0.85,0.15+0.04*(SM2-SM1+1));
1548 hAvECluster->SetAxisRange(hAvECluster->GetMinimum()*0.5, hAvECluster->GetMaximum()*1.1,"Y");
1549 hAvECluster->SetTitleOffset(1.7,"Y");
1550 hAvECluster->SetLineWidth(2);
1551 hAvECluster->Draw();
1552 leg->AddEntry(hAvECluster, "All SM","l");
1553 for (Int_t sm = SM1; sm <= SM2; sm++) {
1554 hAvEClusterSM[sm]->SetLineColor(colorsSM[sm]);
1555 hAvEClusterSM[sm]->SetLineWidth(1);
1556 hAvEClusterSM[sm]->Draw("same");
1557 leg->AddEntry(hAvEClusterSM[sm], Form("SM%i",sm),"l");
1559 hAvECluster->Draw("same"); // to top
1562 // average number of clusters
1564 gPad->SetLeftMargin(0.14);
1565 gPad->SetRightMargin(0.06);
1566 leg = new TLegend(0.65,0.15,0.85,0.15+0.04*(SM2-SM1+1));
1568 hAvNCluster->SetAxisRange(0, hAvNCluster->GetMaximum()*1.3,"Y");
1569 hAvNCluster->SetTitleOffset(1.7,"Y");
1570 hAvNCluster->SetLineWidth(2);
1571 hAvNCluster->Draw();
1572 leg->AddEntry(hAvNCluster, Form("All SM/%i",SM2-SM1+1),"l");
1573 for (Int_t sm = SM1; sm <= SM2; sm++) {
1574 hAvNClusterSM[sm]->SetLineColor(colorsSM[sm]);
1575 hAvNClusterSM[sm]->SetLineWidth(1);
1576 hAvNClusterSM[sm]->Draw("same");
1577 leg->AddEntry(hAvNClusterSM[sm], Form("SM%i",sm),"l");
1579 hAvNCluster->Draw("same"); // to top
1582 // average number of cells in cluster
1584 gPad->SetLeftMargin(0.14);
1585 gPad->SetRightMargin(0.06);
1586 leg = new TLegend(0.65,0.15,0.85,0.15+0.04*(SM2-SM1+1));
1588 hAvNCellsInCluster->SetAxisRange(0, hAvNCellsInCluster->GetMaximum()*1.3,"Y");
1589 hAvNCellsInCluster->SetTitleOffset(1.7,"Y");
1590 hAvNCellsInCluster->SetLineWidth(2);
1591 hAvNCellsInCluster->Draw();
1592 leg->AddEntry(hAvNCellsInCluster, "All SM","l");
1593 for (Int_t sm = SM1; sm <= SM2; sm++) {
1594 hAvNCellsInClusterSM[sm]->SetLineColor(colorsSM[sm]);
1595 hAvNCellsInClusterSM[sm]->SetLineWidth(1);
1596 hAvNCellsInClusterSM[sm]->Draw("same");
1597 leg->AddEntry(hAvNCellsInClusterSM[sm], Form("SM%i",sm),"l");
1599 hAvNCellsInCluster->Draw("same"); // to top
1603 /* Draw the same vs number of events
1606 TCanvas *c2 = new TCanvas(Form("ClusterAveragesEvents%s",suff),
1607 Form("ClusterAveragesEvents%s",suff), 999,333);
1610 // average cluster energy
1612 gPad->SetLeftMargin(0.14);
1613 gPad->SetRightMargin(0.08);
1614 leg = new TLegend(0.65,0.15,0.85,0.15+0.04*(SM2-SM1+1));
1616 TH1* hAvEClusterEvents = ConvertMapRuns2Events(nruns, runNumbers, hAvECluster);
1617 hAvEClusterEvents->Draw();
1618 leg->AddEntry(hAvEClusterEvents, "All SM","l");
1619 for (Int_t sm = SM1; sm <= SM2; sm++) {
1620 TH1* hAvEClusterSMEvents = ConvertMapRuns2Events(nruns, runNumbers, hAvEClusterSM[sm]);
1621 hAvEClusterSMEvents->Draw("same");
1622 leg->AddEntry(hAvEClusterSMEvents, Form("SM%i",sm),"l");
1624 hAvEClusterEvents->Draw("same"); // to top
1627 // average number of clusters
1629 gPad->SetLeftMargin(0.14);
1630 gPad->SetRightMargin(0.08);
1631 leg = new TLegend(0.65,0.15,0.85,0.15+0.04*(SM2-SM1+1));
1633 TH1* hAvNClusterEvents = ConvertMapRuns2Events(nruns, runNumbers, hAvNCluster);
1634 hAvNClusterEvents->Draw();
1635 leg->AddEntry(hAvNClusterEvents, Form("All SM/%i",SM2-SM1+1),"l");
1636 for (Int_t sm = SM1; sm <= SM2; sm++) {
1637 TH1* hAvNClusterSMEvents = ConvertMapRuns2Events(nruns, runNumbers, hAvNClusterSM[sm]);
1638 hAvNClusterSMEvents->Draw("same");
1639 leg->AddEntry(hAvNClusterSMEvents, Form("SM%i",sm),"l");
1641 hAvNClusterEvents->Draw("same"); // to top
1644 // average number of cells in cluster
1646 gPad->SetLeftMargin(0.14);
1647 gPad->SetRightMargin(0.08);
1648 leg = new TLegend(0.65,0.15,0.85,0.15+0.04*(SM2-SM1+1));
1650 TH1* hAvNCellsInClusterEvents = ConvertMapRuns2Events(nruns, runNumbers, hAvNCellsInCluster);
1651 hAvNCellsInClusterEvents->Draw();
1652 leg->AddEntry(hAvNCellsInClusterEvents, "All SM","l");
1653 for (Int_t sm = SM1; sm <= SM2; sm++) {
1654 TH1* hAvNCellsInClusterSMEvents = ConvertMapRuns2Events(nruns, runNumbers, hAvNCellsInClusterSM[sm]);
1655 hAvNCellsInClusterSMEvents->Draw("same");
1656 leg->AddEntry(hAvNCellsInClusterSMEvents, Form("SM%i",sm),"l");
1658 hAvNCellsInClusterEvents->Draw("same"); // to top
1666 //_________________________________________________________________________
1667 Double_t GetAcceptance(Int_t sm, TH2* hmap, Int_t ri)
1669 // Returns [#cells - #dead]/#cells for a supermodule.
1670 // hmap -- dead/noisy cell map;
1673 // guess number of cells per supermodule
1674 Int_t nSM = 1152; // EMCAL
1675 if (hmap->GetXaxis()->GetXmin() == 1) {// PHOS
1677 sm--; // starts from 1, convenient from 0
1682 for (Int_t k = 1; k <= nSM; k++)
1683 if (hmap->GetBinContent(nSM*sm + k, ri+1) < 0)
1686 return ((Double_t) (nSM - ndead))/nSM;
1689 //_________________________________________________________________________
1690 void DrawPi0Slice(Int_t run, Int_t sm = -1)
1692 // Draw the pi0 peak;
1693 // run,sm -- run number and supermodule to take;
1694 // sm < 0 -- draw for whole detector.
1697 if (sm >= 0) {//particular supermodule
1698 h = (TH1*) gFile->Get(Form("run%i_hPi0MassSM%iSM%i",run,sm,sm));
1699 h->SetName(Form("hPi0SliceSM%i_run%i",sm,run));
1700 h->SetTitle(Form("#pi^{0} in SM%i, run %i, %.2fM events", sm, run, GetNumberOfEvents(run)/1e+6));
1702 else {// whole detector
1703 for (Int_t sm1 = 0; sm1 < 10; sm1++)
1704 for (Int_t sm2 = sm1; sm2 < 10; sm2++) {
1705 TH1* one = (TH1*) gFile->Get(Form("run%i_hPi0MassSM%iSM%i",run,sm1,sm2));
1714 h->SetName(Form("hPi0Slice_run%i",run));
1715 h->SetTitle(Form("#pi^{0} in all SM, run %i, %.2fM events", run, GetNumberOfEvents(run)/1e+6));
1718 h->SetXTitle("M_{#gamma#gamma}, GeV");
1719 h->SetYTitle("Entries");
1720 h->SetTitleOffset(1.7,"Y");
1723 if (sm >= 0) c1 = new TCanvas(Form("hPi0SliceSM%i_run%i",sm,run),Form("hPi0SliceSM%i_run%i",sm,run), 400,400);
1724 else c1 = new TCanvas(Form("hPi0Slice_run%i",run),Form("hPi0Slice_run%i",run), 400,400);
1726 gPad->SetLeftMargin(0.14);
1727 gPad->SetRightMargin(0.06);
1730 Double_t nraw, enraw, mass, emass, sigma, esigma;
1731 FitPi0(h, nraw, enraw, mass, emass, sigma, esigma);
1734 TF1* fitfun = h->GetFunction("fitfun");
1736 Double_t emin, emax;
1737 fitfun->GetRange(emin, emax);
1739 backgr = new TF1("mypol2", "[0] + [1]*(x-0.135) + [2]*(x-0.135)^2", emin, emax);
1740 backgr->SetLineColor(kBlue);
1741 backgr->SetLineWidth(2);
1742 backgr->SetLineStyle(3);
1743 backgr->SetParameters(fitfun->GetParameter(3), fitfun->GetParameter(4), fitfun->GetParameter(5));
1744 backgr->Draw("same");
1750 //_________________________________________________________________________
1751 void DrawPi0Averages(Int_t nruns, Int_t runNumbers[], Bool_t samesm = kFALSE, TH2* hmap = NULL)
1753 // Draw average number of pi0s per event, pi0 mass position and pi0 width per run index.
1754 // Errors are also drawn.
1756 // samesm -- take only pi0s for which gammas were is same SM;
1757 // hmap -- if not NULL, do simple (area law based) acceptance correction.
1759 // TODO: PHOS needs pi0 between SMs rather than in one SM
1762 char *suff = TString(Form("%s%s", samesm ? "_sameSM":"", hmap ? "_corr4accep":"")).Data();
1764 // supermodule region
1767 while (SM1 <= SM2 && !gFile->Get(Form("run%i_hPi0MassSM%iSM%i",runNumbers[0],SM1,SM1)) ) SM1++;
1768 while (SM2 >= SM1 && !gFile->Get(Form("run%i_hPi0MassSM%iSM%i",runNumbers[0],SM2,SM2)) ) SM2--;
1770 // initialize histograms for the entire detector
1771 TH1* hPi0Num = new TH1F(Form("hPi0Num%s",suff),"Average number of #pi^{0}s per event", nruns,0,nruns);
1772 hPi0Num->SetXTitle("Run index");
1773 hPi0Num->SetYTitle("Number of #pi^{0}s");
1775 TH1* hPi0Mass = new TH1F(Form("hPi0Mass%s",suff),"#pi^{0} mass position", nruns,0,nruns);
1776 hPi0Mass->SetXTitle("Run index");
1777 hPi0Mass->SetYTitle("M_{#pi^{0}}, GeV");
1779 TH1* hPi0Sigma = new TH1F(Form("hPi0Sigma%s",suff),"#pi^{0} width", nruns,0,nruns);
1780 hPi0Sigma->SetXTitle("Run index");
1781 hPi0Sigma->SetYTitle("#sigma_{#pi^{0}}, GeV");
1783 // initialize histograms per SM
1785 TH1* hPi0MassSM[10];
1786 TH1* hPi0SigmaSM[10];
1788 for (Int_t sm = SM1; sm <= SM2; sm++) {
1789 hPi0NumSM[sm] = new TH1F(Form("hPi0NumSM%i%s",sm,suff),"", nruns,0,nruns);
1790 hPi0MassSM[sm] = new TH1F(Form("hPi0MassSM%i%s",sm,suff),"", nruns,0,nruns);
1791 hPi0SigmaSM[sm] = new TH1F(Form("hPi0SigmaSM%i%s",sm,suff),"", nruns,0,nruns);
1795 for (Int_t ri = 0; ri < nruns; ri++)
1797 fprintf(stderr, "\rDrawPi0Averages(): analysing run index %i/%i ... ", ri, nruns-1);
1799 Int_t nevents = GetNumberOfEvents(runNumbers[ri]);
1802 for (Int_t sm = SM1; sm <= SM2; sm++) {
1803 TH1* h = (TH1*) gFile->Get(Form("run%i_hPi0MassSM%iSM%i",runNumbers[ri],sm,sm));
1805 // supermodule acceptance
1806 Double_t accep = 1.;
1807 if (hmap) accep = GetAcceptance(sm, hmap, ri);
1808 if (accep == 0) continue; // missing SM
1810 Double_t nraw, enraw, mass, emass, sigma, esigma;
1811 FitPi0(h, nraw, enraw, mass, emass, sigma, esigma);
1813 hPi0NumSM[sm]->SetBinContent(ri+1, nraw/nevents/accep);
1814 hPi0MassSM[sm]->SetBinContent(ri+1, mass);
1815 hPi0SigmaSM[sm]->SetBinContent(ri+1, sigma);
1817 hPi0NumSM[sm]->SetBinError(ri+1, enraw/nevents/accep);
1818 hPi0MassSM[sm]->SetBinError(ri+1, emass);
1819 hPi0SigmaSM[sm]->SetBinError(ri+1, esigma);
1822 } // supermodule loop
1825 /* fill for the entire detector
1827 TH1* hsum = (TH1*) gFile->Get(Form("run%i_hPi0MassSM%iSM%i",runNumbers[ri],SM1,SM1));
1828 hsum->SetName("hSumTMP");
1830 // for acceptance correction
1833 for (Int_t sm1 = SM1; sm1 <= SM2; sm1++)
1834 for (Int_t sm2 = sm1; sm2 <= SM2; sm2++) {
1835 if (sm1 == SM1 && sm2 == SM2) continue;
1836 if (samesm && sm1 != sm2) continue;
1838 TH1* h = (TH1*) gFile->Get(Form("run%i_hPi0MassSM%iSM%i",runNumbers[ri],sm1,sm2));
1840 // correct for SM acceptance
1842 Double_t accep = GetAcceptance(sm1, hmap, ri);
1843 if (accep > 0) h->Scale(1./accep);
1851 Double_t nraw, enraw, mass, emass, sigma, esigma;
1852 FitPi0(hsum, nraw, enraw, mass, emass, sigma, esigma);
1854 hPi0Num->SetBinContent(ri+1, nraw/nevents/(SM2-SM1+1-noSM));
1855 hPi0Mass->SetBinContent(ri+1, mass);
1856 hPi0Sigma->SetBinContent(ri+1, sigma);
1858 hPi0Num->SetBinError(ri+1, enraw/nevents/(SM2-SM1+1-noSM));
1859 hPi0Mass->SetBinError(ri+1, emass);
1860 hPi0Sigma->SetBinError(ri+1, esigma);
1865 fprintf(stderr, "\n");
1871 // number of pi0s vs run index
1872 TCanvas *c1 = new TCanvas(Form("hPi0NumRuns%s",suff),Form("hPi0NumRuns%s",suff), 800,400);
1873 gPad->SetLeftMargin(0.08);
1874 gPad->SetRightMargin(0.06);
1875 gPad->SetTopMargin(0.12);
1876 TLegend *leg = new TLegend(0.75,0.15,0.95,0.15+0.04*(SM2-SM1+1));
1878 hPi0Num->SetAxisRange(0., hPi0Num->GetMaximum()*1.2, "Y");
1879 hPi0Num->SetTitleOffset(1.7, "Y");
1880 hPi0Num->SetLineWidth(1);
1882 leg->AddEntry(hPi0Num, Form("All SM/%i",SM2-SM1+1),"l");
1883 for (Int_t sm = SM1; sm <= SM2; sm++) {
1884 hPi0NumSM[sm]->SetLineColor(colorsSM[sm]);
1885 hPi0NumSM[sm]->SetLineWidth(1);
1886 hPi0NumSM[sm]->Draw("same");
1887 leg->AddEntry(hPi0NumSM[sm], Form("SM%i",sm),"l");
1889 hPi0Num->Draw("same"); // to the top
1890 hPi0Num->Draw("hist,same");
1894 // number of pi0s vs event count
1895 TCanvas *c2 = new TCanvas(Form("hPi0NumEvents%s",suff),Form("hPi0NumEvents%s",suff), 800,400);
1896 gPad->SetLeftMargin(0.08);
1897 gPad->SetRightMargin(0.06);
1898 gPad->SetTopMargin(0.12);
1899 leg = new TLegend(0.75,0.15,0.92,0.15+0.04*(SM2-SM1+1));
1901 TH1* hPi0NumEvents = ConvertMapRuns2Events(nruns, runNumbers, hPi0Num);
1902 hPi0NumEvents->Draw();
1903 leg->AddEntry(hPi0NumEvents, Form("All SM/%i",SM2-SM1+1),"l");
1904 for (Int_t sm = SM1; sm <= SM2; sm++) {
1905 TH1* hPi0NumSMEvents = ConvertMapRuns2Events(nruns, runNumbers, hPi0NumSM[sm]);
1906 hPi0NumSMEvents->Draw("same");
1907 leg->AddEntry(hPi0NumSMEvents, Form("SM%i",sm),"l");
1909 hPi0NumEvents->Draw("same"); // to the top
1910 hPi0NumEvents->Draw("hist,same");
1914 // pi0 mass and width vs run index
1915 TCanvas *c3 = new TCanvas(Form("hPi0MassSigmaRuns%s",suff),Form("hPi0MassSigmaRuns%s",suff), 800,400);
1919 gPad->SetLeftMargin(0.16);
1920 gPad->SetRightMargin(0.04);
1921 leg = new TLegend(0.75,0.15,0.95,0.15+0.04*(SM2-SM1+1));
1923 hPi0Mass->SetAxisRange(0.125, 0.145, "Y");
1924 hPi0Mass->SetTitleOffset(2.0, "Y");
1925 hPi0Mass->SetLineWidth(1);
1927 leg->AddEntry(hPi0Mass, "All SM","l");
1928 for (Int_t sm = SM1; sm <= SM2; sm++) {
1929 hPi0MassSM[sm]->SetLineColor(colorsSM[sm]);
1930 hPi0MassSM[sm]->SetLineWidth(1);
1931 hPi0MassSM[sm]->Draw("same");
1932 leg->AddEntry(hPi0MassSM[sm], Form("SM%i",sm),"l");
1934 hPi0Mass->Draw("same"); // to the top
1935 hPi0Mass->Draw("hist,same");
1939 gPad->SetLeftMargin(0.16);
1940 gPad->SetRightMargin(0.04);
1941 leg = new TLegend(0.75,0.15,0.95,0.15+0.04*(SM2-SM1+1));
1943 hPi0Sigma->SetAxisRange(0., hPi0Sigma->GetMaximum()*1.5, "Y");
1944 hPi0Sigma->SetTitleOffset(2.0, "Y");
1945 hPi0Sigma->SetLineWidth(1);
1947 leg->AddEntry(hPi0Sigma, "All SM","l");
1948 for (Int_t sm = SM1; sm <= SM2; sm++) {
1949 hPi0SigmaSM[sm]->SetLineColor(colorsSM[sm]);
1950 hPi0SigmaSM[sm]->SetLineWidth(1);
1951 hPi0SigmaSM[sm]->Draw("same");
1952 leg->AddEntry(hPi0SigmaSM[sm], Form("SM%i",sm),"l");
1954 hPi0Sigma->Draw("same"); // to the top
1955 hPi0Sigma->Draw("hist,same");
1959 // pi0 mass and width vs number of events
1960 TCanvas *c4 = new TCanvas(Form("hPi0MassSigmaEvents%s",suff),Form("hPi0MassSigmaEvents%s",suff), 800,400);
1964 gPad->SetLeftMargin(0.16);
1965 gPad->SetRightMargin(0.08);
1966 leg = new TLegend(0.75,0.15,0.91,0.15+0.04*(SM2-SM1+1));
1968 TH1* hPi0MassEvents = ConvertMapRuns2Events(nruns, runNumbers, hPi0Mass);
1969 hPi0MassEvents->Draw();
1970 leg->AddEntry(hPi0MassEvents, "All SM","l");
1971 for (Int_t sm = SM1; sm <= SM2; sm++) {
1972 TH1* hPi0MassSMEvents = ConvertMapRuns2Events(nruns, runNumbers, hPi0MassSM[sm]);
1973 hPi0MassSMEvents->Draw("same");
1974 leg->AddEntry(hPi0MassSMEvents, Form("SM%i",sm),"l");
1976 hPi0MassEvents->Draw("same"); // to the top
1977 hPi0MassEvents->Draw("hist,same");
1981 gPad->SetLeftMargin(0.16);
1982 gPad->SetRightMargin(0.08);
1983 leg = new TLegend(0.75,0.15,0.91,0.15+0.04*(SM2-SM1+1));
1985 TH1* hPi0SigmaEvents = ConvertMapRuns2Events(nruns, runNumbers, hPi0Sigma);
1986 hPi0SigmaEvents->Draw();
1987 leg->AddEntry(hPi0SigmaEvents, "All SM","l");
1988 for (Int_t sm = SM1; sm <= SM2; sm++) {
1989 TH1* hPi0SigmaSMEvents = ConvertMapRuns2Events(nruns, runNumbers, hPi0SigmaSM[sm]);
1990 hPi0SigmaSMEvents->Draw("same");
1991 leg->AddEntry(hPi0SigmaSMEvents, Form("SM%i",sm),"l");
1993 hPi0SigmaEvents->Draw("same"); // to the top
1994 hPi0SigmaEvents->Draw("hist,same");
2004 //_________________________________________________________________________
2005 void DrawPi0Distributions(char *suff, Int_t nbins = 100)
2007 // Draw distributions for
2008 // 1) average number of pi0s per event;
2009 // 2) pi0 mass position;
2012 // Must be called after DrawPi0Averages() because it
2013 // searches for the corresponding histograms by name.
2015 // suff -- histograms suffix;
2016 // nbins -- number of bins in distributions.
2018 TCanvas *c1 = new TCanvas(Form("Pi0Distributions%s",suff),Form("Pi0Distributions%s",suff), 999,333);
2022 c1->cd(1)->SetLogy();
2023 gPad->SetLeftMargin(0.16);
2024 gPad->SetRightMargin(0.04);
2025 TH1* hPi0Num = (TH1*) gROOT->FindObject(Form("hPi0Num%s",suff));
2026 MakeDistribution(hPi0Num,nbins)->Draw();
2029 c1->cd(2)->SetLogy();
2030 gPad->SetLeftMargin(0.16);
2031 gPad->SetRightMargin(0.04);
2032 TH1* hPi0Mass = (TH1*) gROOT->FindObject(Form("hPi0Mass%s",suff));
2033 MakeDistribution(hPi0Mass,nbins)->Draw();
2036 c1->cd(3)->SetLogy();
2037 gPad->SetLeftMargin(0.16);
2038 gPad->SetRightMargin(0.04);
2039 TH1* hPi0Sigma = (TH1*) gROOT->FindObject(Form("hPi0Sigma%s",suff));
2040 MakeDistribution(hPi0Sigma,nbins)->Draw();
2045 //_________________________________________________________________________
2046 TH1* MakeDistribution(TH1* inhisto, Int_t nbins)
2048 // Returns distribution for the input histogram.
2050 // inhisto -- input histogram;
2051 // nbins -- number of bins in distribution.
2053 // initialize distribution; 1.01 -- to see the last bin
2054 TH1* distr = new TH1F(Form("%sDistr",inhisto->GetName()),"",
2055 nbins, inhisto->GetMinimum(),inhisto->GetMaximum()*1.01);
2056 distr->SetXTitle(inhisto->GetYaxis()->GetTitle());
2057 distr->SetYTitle("Number of runs");
2059 // fill distribution
2060 for (Int_t i = 1; i <= inhisto->GetNbinsX(); i++)
2061 distr->Fill(inhisto->GetBinContent(i));
2066 //_________________________________________________________________________
2067 void FitPi0(TH1* h, Double_t &nraw, Double_t &enraw,
2068 Double_t &mass, Double_t &emass,
2069 Double_t &sigma, Double_t &esigma,
2070 Double_t emin = 0.05, Double_t emax = 0.3, Int_t rebin = 1)
2072 // Fits the pi0 peak with crystal ball + pol2,
2073 // fills number of pi0s, mass, width and their errors.
2079 if (h->GetEntries() == 0) return NULL;
2081 if (rebin > 1) h->Rebin(rebin);
2083 // crystal ball parameters (fixed by hand for EMCAL); TODO: PHOS parameters?
2084 Double_t alpha = 1.1; // alpha >= 0
2085 Double_t n = 2.; // n > 1
2087 // CB tail parameters
2088 Double_t a = pow((n/alpha), n) * TMath::Exp(-alpha*alpha/2.);
2089 Double_t b = n/alpha - alpha;
2091 // integral value under crystal ball with amplitude = 1, sigma = 1
2092 // (will be sqrt(2pi) at alpha = infinity)
2093 Double_t nraw11 = a * pow(b+alpha, 1.-n)/(n-1.) + TMath::Sqrt(TMath::Pi()/2.) * TMath::Erfc(-alpha/TMath::Sqrt(2.));
2095 // signal (crystal ball);
2096 new TF1("cball", Form("(x-[1])/[2] > -%f ? \
2097 [0]*exp(-(x-[1])*(x-[1])/(2*[2]*[2])) \
2098 : [0]*%f*(%f-(x-[1])/[2])^(-%f)", alpha, a, b, n));
2101 new TF1("mypol2", "[0] + [1]*(x-0.135) + [2]*(x-0.135)^2", emin, emax);
2103 // signal + background
2104 TF1 *fitfun = new TF1("fitfun", "cball + mypol2", emin, emax);
2105 fitfun->SetParNames("A", "M", "#sigma", "a_{0}", "a_{1}", "a_{2}");
2106 fitfun->SetLineColor(kRed);
2107 fitfun->SetLineWidth(2);
2109 // make a preliminary fit to estimate parameters
2110 TF1* ff = new TF1("fastfit", "gaus(0) + [3]");
2111 ff->SetParLimits(0, 5., h->GetMaximum()*1.5);
2112 ff->SetParLimits(1, 0.1, 0.2);
2113 ff->SetParLimits(2, 0.005,0.030);
2114 ff->SetParameters(h->GetMaximum()/3., 0.135, 0.010, 0.);
2115 h->Fit(ff, "0QEM", "", 0.105, 0.165);
2117 fitfun->SetParameters(ff->GetParameter(0), ff->GetParameter(1), ff->GetParameter(2), ff->GetParameter(3));
2118 h->Fit(fitfun,"QLEMR", "");
2120 // calculate pi0 values
2121 mass = fitfun->GetParameter(1);
2122 emass = fitfun->GetParError(1);
2124 sigma = fitfun->GetParameter(2);
2125 esigma = fitfun->GetParError(2);
2127 Double_t A = fitfun->GetParameter(0);
2128 Double_t eA = fitfun->GetParError(0);
2130 nraw = nraw11 * A * sigma / h->GetBinWidth(1);
2131 enraw = nraw * (eA/A + esigma/sigma);
2135 /* COMMON FUNCTIONS FOR RUN NUMBERS EXTRACTION, VISUALIZATION AND FILTERING.
2137 * Input: histogram hNEventsProcessedPerRun which is taken from gFile.
2140 //_________________________________________________________________________
2141 void GetRunNumbers(Int_t &nruns, Int_t runNumbers[])
2143 // Fill runNumbers array with run numbers according to hNEventsProcessedPerRun
2144 // histogram: actually analysed runs have positive bin content.
2146 TH1* hNEventsProcessedPerRun = (TH1*) gFile->Get("hNEventsProcessedPerRun");
2148 // check underflow/overflow
2149 if (hNEventsProcessedPerRun->GetBinContent(0) > 0.5)
2150 Warning("GetRunNumbers", "some run numbers are in underflow; they will be absent in the list of runs");
2151 if (hNEventsProcessedPerRun->GetBinContent(hNEventsProcessedPerRun->GetNbinsX()+1) > 0.5)
2152 Warning("GetRunNumbers", "some run numbers are in overflow; they will be absent in the list of runs");
2156 for (Int_t i = 1; i <= hNEventsProcessedPerRun->GetNbinsX(); i++)
2157 if (hNEventsProcessedPerRun->GetBinContent(i) > 0.5) {
2158 runNumbers[nruns] = hNEventsProcessedPerRun->GetBinLowEdge(i);
2163 //_________________________________________________________________________
2164 Int_t GetNumberOfEvents(Int_t run)
2166 // Return number of events in run;
2167 // run -- run number.
2169 TH1* hNEventsProcessedPerRun = (TH1*) gFile->Get("hNEventsProcessedPerRun");
2171 // round the number of events to avoid precision surprizes
2172 return TMath::Nint( hNEventsProcessedPerRun->GetBinContent(hNEventsProcessedPerRun->FindBin(run)) );
2175 //_________________________________________________________________________
2176 Long64_t GetTotalNumberOfEvents(Int_t nruns, Int_t runNumbers[])
2178 // Return total number of events in all runs
2180 Long64_t ntotal = 0;
2182 for (Int_t i = 0; i < nruns; i++)
2183 ntotal += GetNumberOfEvents(runNumbers[i]);
2188 //_________________________________________________________________________
2189 void DrawRunsDistribution(Int_t nruns, Int_t runNumbers[], Int_t dnbins = 100)
2191 // Draw events and events distribution;
2192 // dnbins -- number of bins in distribution.
2194 // initialize run index histogram ...
2195 TH1* hNEventsPerRunIndex = new TH1F("hNEventsPerRunIndex", "Number of processed events per run index", nruns,0,nruns);
2196 hNEventsPerRunIndex->SetXTitle("Run index");
2197 hNEventsPerRunIndex->SetYTitle("Number of events");
2200 for (Int_t i = 0; i < nruns; i++)
2201 hNEventsPerRunIndex->SetBinContent(i+1, GetNumberOfEvents(runNumbers[i]));
2203 // initialize distribution; 1.01 - to see the last bin
2204 TH1* distrib = new TH1F("hNEventsPerRunIndexDistr", "", dnbins, 0, hNEventsPerRunIndex->GetMaximum()*1.01);
2205 distrib->SetXTitle("Number of events");
2206 distrib->SetYTitle("Number of runs");
2208 // fill distribution
2209 for (Int_t i = 1; i <= nruns; i++)
2210 distrib->Fill(hNEventsPerRunIndex->GetBinContent(i));
2212 // draw histogram + distribution
2213 TCanvas *c1 = new TCanvas("hNEventsPerRunIndex","hNEventsPerRunIndex", 800,400);
2217 gPad->SetLeftMargin(0.12);
2218 gPad->SetRightMargin(0.08);
2219 gPad->SetTopMargin(0.12);
2220 hNEventsPerRunIndex->SetTitleOffset(1.7,"Y");
2221 hNEventsPerRunIndex->Draw();
2224 gPad->SetLeftMargin(0.12);
2225 gPad->SetRightMargin(0.08);
2231 //_________________________________________________________________________
2232 void ExcludeSmallRuns(Int_t &nruns, Int_t runNumbers[], Int_t nEventsMin = 100e+3)
2234 // Exclude runs with number of events < nEventsMin
2236 printf("Excluding runs (< %.0fk events):", nEventsMin/1000.);
2238 for (Int_t i = 0; i < nruns; i++)
2239 if (GetNumberOfEvents(runNumbers[i]) < nEventsMin) {
2240 printf(" %i", runNumbers[i]);
2242 runNumbers[i] = runNumbers[nruns];
2248 SortArray(nruns, runNumbers);
2251 //_________________________________________________________________________
2252 void ExcludeRunNumbers(Int_t &nruns, Int_t runNumbers[], Int_t nexcl, Int_t runs2Exclude[])
2254 // Exclude particular runs.
2256 // nexcl -- number of runs in runs2Exclude[];
2257 // runs2Exclude -- array with run numbers to exclude.
2259 for (Int_t i = 0; i < nruns; i++)
2260 for (Int_t e = 0; e < nexcl; e++)
2261 if (runNumbers[i] == runs2Exclude[e]) {
2262 Printf("Excluding run: %i", runs2Exclude[e]);
2264 runNumbers[i] = runNumbers[nruns];
2269 SortArray(nruns, runNumbers);
2272 //_________________________________________________________________________
2273 void SortArray(const Int_t nruns, Int_t runNumbers[])
2276 // used for runNumbers array sorting.
2278 Int_t indices[nruns];
2279 Int_t runNumbers_unsort[nruns];
2281 for (Int_t i = 0; i < nruns; i++)
2282 runNumbers_unsort[i] = runNumbers[i];
2284 TMath::Sort(nruns, runNumbers_unsort, indices, kFALSE);
2286 for (Int_t i = 0; i < nruns; i++)
2287 runNumbers[i] = runNumbers_unsort[indices[i]];
2290 void PrintRunNumbers(Int_t nruns, Int_t runNumbers[])
2292 // Print a table run index / run number / number of events.
2295 Printf("| index | run number | nevents |");
2296 Printf("|-------|------------|-----------|");
2298 for (Int_t i = 0; i < nruns; i++)
2299 Printf("| %-4i | %-9i | %-8i |", i, runNumbers[i], GetNumberOfEvents(runNumbers[i]));
2301 Printf("| Events in total: %-10lli |\n", GetTotalNumberOfEvents(nruns, runNumbers));