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 = "LHC11e_cpass1_CellQA_AnyInt.root";
51 char *infile = "LHC11e_cpass1_CellQA_PHI7.root";
53 Bool_t trendPlotEvents = kFALSE;
56 Color_t colorsSM[] = {0,2,3,4};
58 // cells to exclude from averages calculations and their number: bad cells can
59 // mess up the results, it is suggested to list (known and newly found with
60 // this code) problematic cell candidates explicitly
62 Int_t excells[] = {1};
65 // ====================== End of parameters section ===========================
67 static TFile *gInputFile;
69 void getCellsRunsQAPHOS()
71 // Entry point for the analysis.
74 gStyle->SetOptStat(0);
75 gStyle->SetOptFit(111);
76 gStyle->SetPalette(1);
77 gStyle->SetFillColor(10);
79 Printf("Input: %s", infile);
80 gInputFile = TFile::Open(infile);
82 // You may wish to extract and draw a particular histogram from infile.
83 // Here are the examples:
85 // TH1* hNEventsProcessedPerRun = (TH1*) gInputFile->Get("hNEventsProcessedPerRun");
86 // hNEventsProcessedPerRun->Draw();
89 // TH1* hNTimesInClusterElow = (TH1*) gInputFile->Get("run146686_hCellLocMaxNTimesInClusterElow");
90 // hNTimesInClusterElow->Draw();
93 // TH1* hETotalClusterEhigh = (TH1*) gInputFile->Get("run146686_hCellLocMaxETotalClusterEhigh");
94 // hETotalClusterEhigh->Draw();
97 // Draw a random cell spectrum;
98 // 0.2-1 GeV -- fit region;
99 // 3 GeV -- the spectrum is drawn up to this value, -1 = no limit;
100 // last argument -- histogram name to process, possible values are:
101 // hCellAmplitude, hCellAmplitudeEHigh, hCellAmplitudeNonLocMax or fhCellAmplitudeEhighNonLocMax.
102 DrawCell(1+gRandom->Integer(10752), 0.25, 1., 4., "hCellAmplitude"); // common cell region for EMCAL2010 and for PHOS
104 // Draw a random cell time spectrum
105 // DrawCellTime(1+gRandom->Integer(10752));
108 /* RUN NUMBER SELECTION SECTION
110 * NOTE: at any point below runs are sorted in chronological order.
113 // array with run numbers and their number
114 Int_t runNumbers[10000];
117 // First, fill run numbers ...
118 GetRunNumbers(nruns, runNumbers);
119 Printf("Total number of runs: %i", nruns);
121 // ... draw events distribution ...
122 // (the last argument is number of bins in this distribution)
123 DrawRunsDistribution(nruns, runNumbers, 100);
125 // ... and exclude runs with number of events < 1k.
126 ExcludeSmallRuns(nruns, runNumbers, 100);
128 // You may wish to exclude particular runs:
129 // Int_t runs2Exclude[] = {111222,333444,555666};
130 // ExcludeRunNumbers(nruns, runNumbers, 3, runs2Exclude);
132 Printf("Number of runs to be analysed: %i", nruns);
134 // Finally, print a nice table with run index / run number / number of events.
135 PrintRunNumbers(nruns, runNumbers);
138 /* PER RUN BAD CELLS SEARCHING CRITERIA
140 * Four primary criteria on a per run basis:
141 * 1 and 2: number of times cell was a local maximum in a cluster at low/high energies;
142 * 3 and 4: total cluster energy for a local maximum cell at low/high energies.
145 // Extract the histograms with dead/noisy cell candidates per run.
146 // For each of the four criteria, for each run and for each cell:
147 // 1) calculate cell factor = [cell value]/[average over cells];
148 // 2) mark cell as dead (candidate) if factor <= factorDead (3rd argument below);
149 // 3) mark cell as noisy (candidate) if factor >= factorNoisy (4th argument below).
151 // Factor thresholds are quite wide by default:
152 // low energy criteria are not very sensitive to them,
153 // while high energy criteria are very sensitive due to limited statistics.
155 // The function below also draws histograms with factor distributions in all the runs
156 // for all the cells. It may help to take the decision about dead/noisy factor thresholds.
157 // The last two arguments -- number of bins and maximum X axis value for such histograms.
159 // Output: hBadCellMap[4] -- bad cell candidates per run for each of the four criteria;
160 // X axis -- cell absId, Y axis -- run index, content: 0=not bad, 1=noisy, -1=dead.
163 TH2** hBadCellMap = FindDeadNoisyCellsPerRun(nruns, runNumbers, 0.05, 2.5, 200, 10.);
165 // Look for noisy channels only:
166 // TH2** hBadCellMap = FindDeadNoisyCellsPerRun(nruns, runNumbers, -1.0, 2.5, 200, 20.);
168 // Look for dead channels only:
169 // TH2** hBadCellMap = FindDeadNoisyCellsPerRun(nruns, runNumbers, 0.05, 1.e9, 200, 10.);
171 // Print cell numbers suggested for exclusion from averages calculation;
172 // see excells[] array in the parameters section in the beginning of the macro.
173 // It is highly suggested to run this macro several times and add the output
174 // of this function to the parameters section.
175 PrintCellsToExcludeFromAverages(hBadCellMap);
177 // The results can be visualized (second argument -- canvas name):
178 // either per run ...
179 DrawDeadNoisyCellMap(hBadCellMap[0], "hMapRuns_NTimesInClusterElow");
180 DrawDeadNoisyCellMap(hBadCellMap[1], "hMapRuns_NTimesInClusterEhigh");
181 // DrawDeadNoisyCellMap(hBadCellMap[2], "hMapRuns_ETotalClusterElow");
182 // DrawDeadNoisyCellMap(hBadCellMap[3], "hMapRuns_ETotalClusterEhigh");
184 // ... or per number of events ...
185 // DrawDeadNoisyCellMap((TH2*)ConvertMapRuns2Events(nruns,runNumbers,hBadCellMap[0]), "hMapEvents_NTimesInClusterElow");
186 // DrawDeadNoisyCellMap((TH2*)ConvertMapRuns2Events(nruns,runNumbers,hBadCellMap[1]), "hMapEvents_NTimesInClusterEhigh");
187 // DrawDeadNoisyCellMap((TH2*)ConvertMapRuns2Events(nruns,runNumbers,hBadCellMap[2]), "hMapEvents_ETotalClusterElow");
188 // DrawDeadNoisyCellMap((TH2*)ConvertMapRuns2Events(nruns,runNumbers,hBadCellMap[3]), "hMapEvents_ETotalClusterEhigh");
190 // ... or we can also draw the percent for each cell being dead/noisy
191 // DrawDeadNoisyCellPercent(nruns, runNumbers, hBadCellMap[0], "hPercent_NTimesInClusterElow");
192 // DrawDeadNoisyCellPercent(nruns, runNumbers, hBadCellMap[1], "hPercent_NTimesInClusterEhigh");
193 // DrawDeadNoisyCellPercent(nruns, runNumbers, hBadCellMap[2], "hPercent_ETotalClusterElow");
194 // DrawDeadNoisyCellPercent(nruns, runNumbers, hBadCellMap[3], "hPercent_ETotalClusterEhigh");
196 // Our main criteria for analytical finding of bad cells is based on the following.
197 // Note that factor distributions for NTimesInCluster and ETotalCluster are very
198 // similar, both at low and at high energy. Note also that we can say the same
199 // for dead/noisy cell visualizations above: they are very similar. This suggests
200 // that NTimesInCluster and ETotalCluster give the same results. The criteria
201 // NTimesInCluster, however, is more calibration-independent (though energy thresholds
202 // are still calibration-dependent) and thus is more reliable and clear. Thus, we
203 // limit ourselves to NTimesInClusterElow/NTimesInClusterEhigh criteria.
204 // Now, you could probably noted that NTimesInClusterEhigh give more dead
205 // cells than that at low energy. This is an expected statistical effect: we have
206 // much smaller number of clusters at high energy. Consequently, we will not use dead
207 // cell candidates at high energy.
209 // Finally, we combine candidates from low/high energies and produce one TH2
210 // histogram which is the primary source of our analytical results.
212 TH2* hBadCellMapPrimary = CombineDeadNoisyElowEhigh(hBadCellMap[0], hBadCellMap[1]);
214 // Note for PHOS: if you are not happy with NTimesInClusterEhigh results (due to a lack of statistics)
215 // uncomment this line:
216 // hBadCellMapPrimary = hBadCellMap[0];
218 // Draw everything for combined
219 DrawDeadNoisyCellMap(hBadCellMapPrimary, "Primary_hMapRuns");
220 // DrawDeadNoisyCellMap((TH2*)ConvertMapRuns2Events(nruns,runNumbers,hBadCellMapPrimary), "Primary_hMapEvents");
221 // DrawDeadNoisyCellPercent(nruns, runNumbers, hBadCellMapPrimary, "Primary_hPercent");
223 // Print full information on cells which are dead/noisy;
224 // kTRUE -- print also the percentage % dead/% noisy
225 PrintDeadNoisyCells(hBadCellMapPrimary, 0.9, 1.); // in 90-100% of runs
226 // PrintDeadNoisyCells(hBadCellMapPrimary, 0.5, 0.9, kTRUE); // in 50-90% of runs
227 // PrintDeadNoisyCells(hBadCellMapPrimary, 0.3, 0.5); // in 30-50% of runs
228 // PrintDeadNoisyCells(hBadCellMapPrimary, 0.1, 0.5); // in 10-50% of runs
230 // visualize dead/noisy cell map for EMCAL/PHOS; requires aliroot
231 DrawOccupancy(nruns, runNumbers, hBadCellMapPrimary, "hDeadNoisyCellsOccupancy");
233 // EMCAL: print full information on missing/noisy parts (e.g. RCUs); requires aliroot
234 // PrintEMCALProblematicBlocks(nruns, runNumbers, hBadCellMapPrimary);
237 /* RUNS QUALITY SECTION: CLUSTER AVERAGES + PI0 AVERAGES
241 // First, draw cluster averages per run;
242 // 1 = minimum number of cells in a cluster;
243 // 0.3GeV = minimum cluster energy;
244 // -1 = maximum cluster energy = infinity (in fact, 20GeV in the default configuration).
245 DrawClusterAveragesPerRun(nruns, runNumbers, 1, 0.3, -1);
247 // Second, draw the same cluster averages per run, but corrected for detector acceptance
248 DrawClusterAveragesPerRun(nruns, runNumbers, 1, 0.3, -1, hBadCellMapPrimary);
250 // Draw random slices of the pi0 peak in one supermodule and in whole detector
251 // Int_t rndRun = runNumbers[gRandom->Integer(nruns)];
252 Int_t rndRun = runNumbers[30];
253 DrawPi0Slice(rndRun, 1);
254 DrawPi0Slice(rndRun, 2);
255 DrawPi0Slice(rndRun, 3);
256 DrawPi0Slice(rndRun, -1);
258 // Draw number of pi0s per event, pi0 mass and width
259 DrawPi0Averages(nruns, runNumbers);
261 // Draw pi0 values per event with pi0 photons both in the same supermodule
262 // DrawPi0Averages(nruns, runNumbers, kTRUE);
264 // Draw pi0 values per event with pi0 photons both in the same supermodule
265 // + correct for supermodule acceptance (should not be treated to be 100% reliable!)
266 // DrawPi0Averages(nruns, runNumbers, kTRUE, hBadCellMap[0]);
268 // Draw pi0 values per event with pi0
269 // + correct for supermodule acceptance (should not be treated to be 100% reliable!)
270 // DrawPi0Averages(nruns, runNumbers, kFALSE, hBadCellMap[0]);
272 // Draw pi0 distributions (helps to take decision on bad runs);
273 // first argument -- suffix from canvas title;
274 // second argument -- number of bins in distributions
275 // DrawPi0Distributions("", 100);
276 // DrawPi0Distributions("SM1", 100)
277 // DrawPi0Distributions("_sameSM", 100);
278 // DrawPi0Distributions("_sameSM_corr4accep", 100);
285 * Lazy/curious boundary is here! -------------------------
286 * Do not pass it if you do not fill curious enough! ;)
288 * The shape analysis below belongs to the main bad cell searching criteria.
289 * However, it may require some extra work in order to make it usable.
291 * The idea behind is simple: fit the energy spectrum of a cell with
292 * the function A*exp(-B*x)/x^2 (which proved to be a very good fit),
293 * draw parameters A, B and chi2/ndf distributions and notice
294 * cells which deviate too much from the averages.
296 * Huge statistics (>20M events) is necessary for this to work reliably.
301 // The analysis below defines a cell as being bad simply if it is outside
302 // of some good region for any of the parameters A, B, chi2/ndf.
303 // The regions are depicted by vertical orange lines. The problem is that these
304 // regions are usually not automatically determined correctly.
306 // The function syntax is the following:
307 // 0.1-1.0 GeV -- fitting region;
308 // hCellAmplitude -- the histogram which to take for processing,
309 // the other possible choice is hCellAmplitudeNonLocMax;
310 // 1000 -- number of bins is distributions.
312 // The next three groups of parameters are:
313 // <text label> <maximum value on distribution> <left edge of the good region> <right edge of the good region>
315 // It is left/right edges which usually require manual tuning.
316 // -1 means to initialize a parameter automatically.
318 TestCellShapes(0.25, 1., "hCellAmplitude", 1000,
319 // maxval / left margin / right margin
322 -1, -1,-1); // fit chi2/ndf
325 /* DISTRIBUTION ANALYSIS
327 * Test for bad cells by plotting a distribution among cells.
329 * It is especially useful in searching for miscalibrated cells.
330 * The function parameters are similar to parameters from shape analysis section.
333 TestCellsTH1(nruns, runNumbers, "hCellLocMaxNTimesInClusterElow",
334 "Number of times cell was local maximum, low energy, per event", "Entries",
335 400, -1,-1,-1); // nbins / maxval / left margin / right margin
337 TestCellsTH1(nruns, runNumbers, "hCellLocMaxNTimesInClusterEhigh",
338 "Number of times cell was local maximum, high energy, per event", "Entries",
339 400, -1,-1,-1); // nbins / maxval / left margin / right margin
341 // TestCellsTH1(nruns, runNumbers, "hCellLocMaxETotalClusterElow",
342 // "Total cluster energy for local maximum cell, low energy, per event", "Entries",
343 // 400, -1,-1,-1); // nbins / maxval / left margin / right margin
345 // TestCellsTH1(nruns, runNumbers, "hCellLocMaxETotalClusterEhigh",
346 // "Total cluster energy for local maximum cell, high energy, per event", "Entries",
347 // 400, -1,-1,-1); // nbins / maxval / left margin / right margin
349 // Three more tests for bad cells:
350 // 1) total deposited energy;
351 // 2) total number of entries;
352 // 3) average energy = [total deposited energy]/[total number of entries].
354 TestCellEandN(0.1, "hCellAmplitude", 1000,
355 // maxval / left margin / right margin
358 -1,-1,-1); // cell E/N
360 // the same at high energies
361 TestCellEandN(0.1, "hCellAmplitudeEhigh", 1000,
362 // maxval / left margin / right margin
365 -1,-1,-1); // cell E/N
368 // The same as above, but for not a local maximum cells; require more statistics
369 // TestCellsTH1(nruns, runNumbers, "hCellNonLocMaxNTimesInClusterElow",
370 // "Number of times cell wasn't local maximum, low energy, per event", "Entries",
371 // 400, -1,-1,-1); // nbins / maxval / left margin / right margin
373 // TestCellsTH1(nruns, runNumbers, "hCellNonLocMaxNTimesInClusterEhigh",
374 // "Number of times cell wasn't local maximum, high energy, per event", "Entries",
375 // 400, -1,-1,-1); // nbins / maxval / left margin / right margin
377 // TestCellsTH1(nruns, runNumbers, "hCellNonLocMaxETotalClusterElow",
378 // "Total cluster energy for not local maximum cell, low energy, per event", "Entries",
379 // 400, -1,-1,-1); // nbins / maxval / left margin / right margin
381 // TestCellsTH1(nruns, runNumbers, "hCellNonLocMaxETotalClusterEhigh",
382 // "Total cluster energy for not local maximum cell, high energy, per event", "Entries",
383 // 400, -1,-1,-1); // nbins / maxval / left margin / right margin
385 // TestCellEandN(0.1, "hCellAmplitudeNonLocMax", 1000,
386 // // maxval / left margin / right margin
387 // -1,-1,-1, // cell E
388 // -1,-1,-1, // cell N
389 // -1,-1,-1); // cell E/N
391 // TestCellEandN(0.1, "hCellAmplitudeEhighNonLocMax", 1000,
392 // // maxval / left margin / right margin
393 // -1,-1,-1, // cell E
394 // -1,-1,-1, // cell N
395 // -1,-1,-1); // cell E/N
405 /* BAD CELLS SEARCHING FUNCTIONS
409 //_________________________________________________________________________
410 TH2** FindDeadNoisyCellsPerRun(const Int_t nruns, Int_t runNumbers[],
411 Double_t factorDead = 0.01, Double_t factorNoisy = 4.,
412 Int_t fnbins = 200, Double_t fxmax = 10.)
414 // Return bad cell candidate maps for the four criteria;
415 // X axis -- cell absId, Y axis -- run index, content: 0=not bad, 1=noisy, -1=dead.
417 // For each run and each cell calculate factor = [cell value]/[averate over cells],
418 // mark cell as noisy if factor >= factorNoisy, mark cell as dead if factor <= factorDead.
420 // fnbins, fxmax -- number of bins and X axis maximum value for factor distributions.
422 // Four criteria in order:
423 char *hname[4] = {"hCellLocMaxNTimesInClusterElow", "hCellLocMaxNTimesInClusterEhigh",
424 "hCellLocMaxETotalClusterElow", "hCellLocMaxETotalClusterEhigh"};
426 TH1* hFactorDistr[4];
427 TH2** hBadCellMap = new TH2*[4];
429 // take one histogram to get binning parameters
430 TH1* one = (TH1*) gInputFile->Get(Form("run%i_%s",runNumbers[0],hname[0]));
431 Int_t ncells = one->GetNbinsX();
432 Double_t amin = one->GetXaxis()->GetXmin();
433 Double_t amax = one->GetXaxis()->GetXmax();
435 // find dead/noisy cell candidates
436 for (Int_t k = 0; k < 4; k++) {
437 hBadCellMap[k] = new TH2C(Form("hBadCellMap_%s_fdead=%.3f_fnoisy=%.1f",hname[k],factorDead,factorNoisy),
438 Form("Dead/noisy cell map (%s, fdead=%.3f, fnoisy=%.1f)",hname[k],factorDead,factorNoisy),
439 ncells,amin,amax, nruns,0,nruns);
440 hBadCellMap[k]->SetXTitle("AbsId");
441 hBadCellMap[k]->SetYTitle("Run index");
442 hBadCellMap[k]->SetTitleOffset(1.7,"Y");
444 hFactorDistr[k] = new TH1F(Form("hFactorDistr_%s_fdead=%.3f_fnoisy=%.1f",
445 hname[k],factorDead,factorNoisy), "", fnbins,0,fxmax);
446 hFactorDistr[k]->SetTitle("Factor distributions in all runs");
447 hFactorDistr[k]->SetXTitle("Factor");
448 hFactorDistr[k]->SetYTitle("Entries");
451 for (Int_t ri = 0; ri < nruns; ri++) {
452 TH1* one = (TH1*) gInputFile->Get(Form("run%i_%s", runNumbers[ri], hname[k]));
454 Error("FindDeadNoisyCellsPerRun","Run %d does contain histogram %s",runNumbers[ri],hname[k]);
458 Double_t av = 0; // average
459 Int_t count = 0; // counted cells
460 for (Int_t c = 1; c <= ncells; c++) {
461 // do not count cells with zero content
462 if (one->GetBinContent(c) == 0) continue;
463 // do not count cells listed in the parameters section in the beginning of the macro
464 if (IsCellMarked4Exclude(one->GetBinLowEdge(c))) continue;
466 av += one->GetBinContent(c);
469 // division by zero checks
471 Warning("FindDeadNoisyCellsPerRun", Form("No cells counted at ri=%i",ri));
477 Warning("FindDeadNoisyCellsPerRun", Form("Average is zero at ri=%i",ri));
481 // find dead/noisy candidates
482 for (Int_t c = 1; c <= ncells; c++) {
483 Double_t fac = one->GetBinContent(c)/av;
484 hFactorDistr[k]->Fill(fac);
486 if (fac <= factorDead)
487 hBadCellMap[k]->SetBinContent(c, ri+1, -1);
488 else if (fac >= factorNoisy)
489 hBadCellMap[k]->SetBinContent(c, ri+1, 1);
497 // draw factor distributions ...
498 TCanvas *c1 = new TCanvas("hFactorDistr", "hFactorDistr", 400,400);
499 gPad->SetLeftMargin(0.12);
500 gPad->SetRightMargin(0.08);
502 hFactorDistr[0]->SetLineColor(kBlue);
503 hFactorDistr[1]->SetLineColor(kRed);
504 hFactorDistr[2]->SetLineColor(kGreen);
505 hFactorDistr[3]->SetLineColor(kOrange);
506 hFactorDistr[0]->Draw();
507 hFactorDistr[1]->Draw("same");
508 hFactorDistr[2]->Draw("same");
509 hFactorDistr[3]->Draw("same");
512 TLegend *leg = new TLegend(0.45,0.65,0.90,0.85);
513 leg->AddEntry(hFactorDistr[0], "NTimesInCluster, low energy","l");
514 leg->AddEntry(hFactorDistr[2], "ETotalCluster, low energy","l");
515 leg->AddEntry(hFactorDistr[1], "NTimesInCluster, high energy","l");
516 leg->AddEntry(hFactorDistr[3], "EtotalCluster, high energy","l");
524 //_________________________________________________________________________
525 void PrintCellsToExcludeFromAverages(TH2** hBadCellMap)
527 // Print cells suggested for exclusion from averages calculation
529 Int_t ncells = hBadCellMap[0]->GetNbinsX();
530 Int_t nruns = hBadCellMap[0]->GetNbinsY();
532 Int_t *suggested = new Int_t[ncells];
533 memset(suggested, 0, ncells*sizeof(Int_t));
535 for (Int_t c = 1; c <= ncells; c++)
536 for (Int_t ri = 1; ri <= nruns; ri++)
537 if (hBadCellMap[0]->GetBinContent(c, ri) != 0 || hBadCellMap[2]->GetBinContent(c, ri) != 0 ||
538 hBadCellMap[1]->GetBinContent(c, ri) > 0 || hBadCellMap[3]->GetBinContent(c, ri) > 0) // NOTE: dead not counted
541 printf("Suggested cells to switch off in averages calculations (copai 'n paste!):\n");
542 printf("Int_t excells[] = {");
544 for (Int_t c = 1; c <= ncells; c++)
545 if (suggested[c-1] >= 0.4*nruns) {// 40% of runs threshold
546 printf("%s%i", n == 0 ? "" : ",", c);
549 printf("};\nInt_t nexc = %i;\n\n",n);
552 //_________________________________________________________________________
553 Bool_t IsCellMarked4Exclude(Int_t absId)
555 // Return true if a cell is in excells[] array
557 static TH1* hExclCells = NULL;
559 // one time initialization
561 hExclCells = new TH1F("hExclCells", "", 20000,0,20000);
562 for (Int_t c = 0; c < nexc; c++)
563 hExclCells->SetBinContent(hExclCells->FindBin(excells[c]), 1);
566 return (hExclCells->GetBinContent(hExclCells->FindBin(absId)) > 0 ? kTRUE : kFALSE);
569 //_________________________________________________________________________
570 void DrawDeadNoisyCellMap(TH2* hmap, char* cname)
572 // Visualize dead/noisy cell map;
573 // cname -- canvas name.
575 // Define only 3 color: dead=blue, noisy=red, normal=white
576 Int_t color_index[]={kBlue,kWhite,kRed};
577 gStyle->SetPalette(3,color_index);
579 TCanvas *c1 = new TCanvas(cname, cname, 0,0,1000,600);
583 // Use hardware numeration of modules (2,3,4)
584 // instead of offline numeration (3,2,1)
585 for (Int_t iM=1; iM<=3; iM++) {
587 gPad->SetLeftMargin(0.04);
588 gPad->SetRightMargin(0.02);
589 gPad->SetTopMargin(0.10);
590 gPad->SetBottomMargin(0.13);
591 hDummy[iM-1] = (TH2*) hmap->Clone(Form("hDummy_%s_mod%d",cname,iM));
592 hDummy[iM-1]->SetTitle(Form("Module %d",5-iM));
593 hDummy[iM-1]->SetLabelSize(0.08,"XY");
594 hDummy[iM-1]->SetTitleSize(0.08,"XY");
595 hDummy[iM-1]->SetTitleOffset(0.8,"X");
596 hDummy[iM-1]->SetTitleOffset(0.25,"Y");
597 hDummy[iM-1]->SetTickLength(0.01,"Y");
598 hDummy[iM-1]->SetNdivisions(520,"X");
599 hDummy[iM-1]->SetAxisRange(3584*(iM-1)+1, 3584*iM, "X");
600 hDummy[iM-1]->Draw("col");
606 // //_________________________________________________________________________
607 // void DrawDeadNoisyCellMap(TH2* hmap, char* cname)
609 // // Visualize dead/noisy cell map;
610 // // cname -- canvas name.
612 // TCanvas *c1 = new TCanvas(cname, cname, 0,0,600,600);
613 // gPad->SetLeftMargin(0.14);
614 // gPad->SetRightMargin(0.06);
616 // // draw dummy to initialize axis ranges
617 // TH2* hDummy = (TH2*) hmap->Clone(Form("hDummy_%s",cname));
621 // for (Int_t c = 1; c <= hmap->GetNbinsX(); c++) //cell number
622 // for (Int_t ri = 1; ri <= hmap->GetNbinsY(); ri++) { //run index
623 // Double_t stat = hmap->GetBinContent(c, ri); // cell status
624 // if (stat == 0) continue;
626 // Double_t x = hmap->GetXaxis()->GetBinCenter(c);
627 // Double_t y1 = hmap->GetYaxis()->GetBinLowEdge(ri);
628 // Double_t y2 = hmap->GetYaxis()->GetBinUpEdge(ri);
630 // // draw a line; FIXME: what is a better choice?
631 // TLine* line = new TLine(x,y1,x,y2);
632 // line->SetLineWidth(1);
633 // if (stat > 0) line->SetLineColor(kRed); // noisy cell
634 // else line->SetLineColor(kBlue); // dead cell
641 //_________________________________________________________________________
642 void DrawDeadNoisyCellPercent(Int_t nruns, Int_t runNumbers[], TH2* hmap, char* cname)
644 // Show percent of runs/events for each cell being problematic;
645 // cname -- canvas name.
647 // binning parameters
648 Int_t ncells = hmap->GetNbinsX();
649 Double_t amin = hmap->GetXaxis()->GetXmin();
650 Double_t amax = hmap->GetXaxis()->GetXmax();
652 // number of times cell was dead/noisy;
653 Int_t *ndeadRuns = new Int_t[ncells];
654 Int_t *nnoisyRuns = new Int_t[ncells];
655 Double_t *ndeadEvents = new Double_t[ncells];
656 Double_t *nnoisyEvents = new Double_t[ncells];
659 for (Int_t c = 0; c < ncells; c++) {
665 for (Int_t ri = 0; ri < nruns; ri++) {
666 Double_t stat = hmap->GetBinContent(c+1, ri+1); // cell status
667 Int_t nevents = GetNumberOfEvents(runNumbers[ri]);
671 nnoisyEvents[c] += nevents;
675 ndeadEvents[c] += nevents;
680 // total number of events
681 Double_t ntotal = GetTotalNumberOfEvents(nruns, runNumbers);
683 TCanvas *c1 = new TCanvas(cname, cname, 0,0,600,600);
684 gPad->SetLeftMargin(0.14);
685 gPad->SetRightMargin(0.06);
687 // draw dummy histogram to initialize canvas
688 TH1* hDummy = new TH1F(Form("hDummy_%s",cname), hmap->GetTitle(), ncells,amin,amax);
689 hDummy->SetAxisRange(0,100, "Y");
690 hDummy->SetXTitle("AbsId");
691 hDummy->SetYTitle("Percent");
700 // draw results, FIXME: is where a better way?
701 for (Int_t c = 0; c < ncells; c++) {
702 Double_t x = hmap->GetXaxis()->GetBinCenter(c+1);
703 Double_t y1 = 100.*ndeadEvents[c]/ntotal;
704 Double_t y2 = 100.*ndeadRuns[c]/nruns;
707 if (ndeadEvents[c] > 0) {
708 line1 = new TLine(x, 0, x, y1);
709 line1->SetLineWidth(1);
710 line1->SetLineColor(kBlue);
715 if (nnoisyEvents[c] > 0) {
716 line2 = new TLine(x, y1, x, y1 + 100.*nnoisyEvents[c]/ntotal);
717 line2->SetLineWidth(1);
718 line2->SetLineColor(kRed);
723 if (ndeadRuns[c] > 0) {
724 line3 = new TLine(x, 0, x, y2);
725 line3->SetLineWidth(1);
726 line3->SetLineStyle(7);
727 line3->SetLineColor(kBlack);
732 if (nnoisyRuns[c] > 0) {
733 line4 = new TLine(x, y2, x, y2 + 100.*nnoisyRuns[c]/nruns);
734 line4->SetLineWidth(1);
735 line4->SetLineStyle(7);
736 line4->SetLineColor(kOrange);
742 TLegend *leg = new TLegend(0.65,0.65,0.92,0.75);
743 if (line1) leg->AddEntry(line1, "% of events, dead","l");
744 if (line2) leg->AddEntry(line2, "% of events, noisy","l");
745 if (line3) leg->AddEntry(line3, "% of runs, dead","l");
746 if (line4) leg->AddEntry(line4, "% of runs, noisy","l");
752 //_________________________________________________________________________
753 TH1* ConvertMapRuns2Events(const Int_t nruns, Int_t runNumbers[], TH1* inhisto)
755 // Returns a histogram in which run index axis is converted to number of events
756 // by making a variable axis bin width.
757 // If inhisto inherits from TH2, convert Y axis; convert X axis otherwise.
760 Double_t *nevents = new Double_t[nruns+1];
762 for (Int_t ri = 0; ri < nruns; ri++)
763 nevents[1+ri] = nevents[ri] + GetNumberOfEvents(runNumbers[ri]);
765 TH1* outh = (TH1*) inhisto->Clone(Form("%s_Events",inhisto->GetName()));
767 if (inhisto->InheritsFrom("TH2")) {
768 outh->GetYaxis()->Set(nruns, nevents);
769 outh->SetYTitle("Number of events");
770 outh->SetTitleOffset(1.7,"Y");
773 outh->GetXaxis()->Set(nruns, nevents);
774 outh->SetXTitle("Number of events");
780 //_________________________________________________________________________
781 TH2* CombineDeadNoisyElowEhigh(TH2* hmapElow, TH2* hmapEhigh)
783 // Combine two maps at low/high energy into one,
784 // do not count dead map at high energy.
785 // NOTE: if cell is dead at Elow and noisy and Ehigh, set content = 2.
787 TH2* hmap_combined = (TH2*) hmapElow->Clone(Form("%s+%s",hmapElow->GetName(),hmapEhigh->GetName()));
789 for (Int_t c = 1; c <= hmapElow->GetNbinsX(); c++)
790 for (Int_t ri = 1; ri <= hmapElow->GetNbinsY(); ri++)
791 if (hmapEhigh->GetBinContent(c, ri) > 0) {
792 hmap_combined->SetBinContent(c, ri, 1);
794 if (hmapElow->GetBinContent(c, ri) < 0)
795 hmap_combined->SetBinContent(c, ri, 2);
798 hmap_combined->SetTitle("Dead/noisy cell map, combined");
800 return hmap_combined;
803 //_________________________________________________________________________
804 void PrintDeadNoisyCells(TH2* hmap, Double_t percent1 = 0.95, Double_t percent2 = 1., Bool_t kPrintPercentage = kFALSE)
806 // Print full information on dead/noisy cells which are present in
807 // percent1-percent2 portion of runs (percent1 excluded, percent2 included).
809 Int_t ncells = hmap->GetNbinsX();
810 Int_t nruns = hmap->GetNbinsY();
814 printf("Dead/noisy cells in >%.1f%% and <=%.1f%% of runs:", 100*percent1, 100*percent2);
816 for (Int_t c = 1; c <= ncells; c++) {
817 Int_t nrdead = 0; // count number of runs for the current cell
820 for (Int_t ri = 1; ri <= nruns; ri++)
821 if (hmap->GetBinContent(c,ri) < 0) nrdead++;
822 else if (hmap->GetBinContent(c,ri) > 0) nrnoisy++;
824 if (nrdead+nrnoisy > percent1*nruns && nrdead+nrnoisy <= percent2*nruns) {
825 printf(" %.0f", hmap->GetBinLowEdge(c));
826 if (kPrintPercentage)
827 printf("(%-2.0f/%-2.0f)", 100.*nrdead/nruns, 100.*nrnoisy/nruns);
832 printf(" (%i total)\n\n", nbad);
835 //_________________________________________________________________________
836 void DrawOccupancy(Int_t nruns, Int_t runNumbers[], TH2* hmap, char* cname)
838 // Draw bad cell map for EMCAL or PHOS;
839 // cname -- canvas name.
841 gStyle->SetPalette(1);
844 if (hmap->GetNbinsX() % 1152 == 0)
845 DrawEMCALOccupancy(nruns, runNumbers, hmap, cname);
847 DrawPHOSOccupancy(nruns, runNumbers, hmap, cname);
850 //_________________________________________________________________________
851 void DrawEMCALOccupancy(Int_t nruns, Int_t runNumbers[], TH2* hmap, char* cname)
853 // Draw bad cell map for EMCAL;
854 // cname -- canvas name.
856 Int_t nmods = hmap->GetNbinsX()/1152; // number of supermodules
857 Int_t vsize = ceil(nmods/2.); // vertical size in canvas
858 Int_t lastSector = (nmods-1)/2; // for pad number calculation
860 TCanvas *c1 = new TCanvas(cname, cname, 800,200*vsize);
861 c1->Divide(2, vsize);
863 for (Int_t sm = 0; sm < nmods; sm++)
865 // top left is SM0, bottom right is SM9
868 c1->cd((lastSector - sector)*2 + side + 1);
870 TH2* hSM = new TH2C(Form("hEMCALSM%i_%s",sm,cname), Form("SM%i",sm), 48,0,48, 24,0,24);
871 hSM->SetXTitle("iEta");
872 hSM->SetYTitle("iPhi");
874 // loop over supermodule cells
875 for (Int_t c = 0; c < 1152; c++) {
876 Int_t absId = 1152 * sm + c;
878 for (Int_t ri = 0; ri < nruns; ri++) {
879 if (hmap->GetBinContent(hmap->FindBin(absId,ri)) == 0) continue;
881 Int_t nModule, eta, phi;
882 AbsId2EtaPhi(absId, nModule, eta, phi);
888 } // supermodule loop
893 //_________________________________________________________________________
894 void DrawPHOSOccupancy(Int_t nruns, Int_t runNumbers[], TH2* hmap, char* cname)
896 // Draw bad cell map for PHOS;
897 // cname -- canvas name.
899 Int_t nmods = hmap->GetNbinsX()/3584; // number of supermodules
900 Int_t vsize = nmods; // vertical size in canvas
902 TCanvas *c1 = new TCanvas(cname, cname, 64*10*vsize,56*10);
905 TFile *fBadMap = TFile::Open("PHOS_BadMap.root","recreate");
906 for (Int_t sm = 1; sm <= nmods; sm++)
909 gPad->SetLeftMargin(0.10);
910 gPad->SetRightMargin(0.15);
911 gPad->SetTopMargin(0.05);
912 gPad->SetBottomMargin(0.10);
914 TH2* hSM = new TH2C(Form("hPHOSSM%i_%s",sm,cname), Form("Module %i",5-sm), 64,0.5,64.5, 56,0.5,56.5);
915 hSM->SetXTitle("x_{cell}");
916 hSM->SetYTitle("z_{cell} ");
917 hSM->SetZTitle("N_{runs} ");
919 // loop over supermodule cells
920 for (Int_t c = 1; c <= 3584; c++) {
921 Int_t absId = 3584*(sm-1) + c;
923 for (Int_t ri = 0; ri < nruns; ri++) {
924 if (hmap->GetBinContent(hmap->FindBin(absId,ri)) == 0) continue;
926 Int_t nModule, xCell, zCell;
927 AbsId2EtaPhi(absId, nModule, xCell, zCell, 1);
928 hSM->Fill(xCell,zCell);
932 hSM->DrawCopy("colz");
933 } // supermodule loop
939 //_________________________________________________________________________
940 void PrintEMCALProblematicBlocks(Int_t nruns, Int_t runNumbers[], TH2* hmap)
942 // Print, on a per run basis, complete information about EMCAL missing
943 // (or dead or noisy) blocks. Missing/noisy EMCAL atomic part is a 2x2
944 // block (288 blocks per supermodule).
946 // number of supermodules
947 Int_t nmods = hmap->GetNbinsX()/1152;
949 Printf("Problematic (missing/dead/noisy) 2x2 block numbers in EMCAL:");
952 for (Int_t ri = 0; ri < nruns; ri++) {
953 Printf("Run %i:", runNumbers[ri]);
956 for (Int_t sm = 0; sm < nmods; sm++) {
958 // will be filled with the number of missing cells (from 0 to 4)
960 memset(blk2x2, 0, 288*sizeof(Int_t));
962 for (Int_t c = 0; c < 1152; c++) {
963 Int_t absId = 1152 * sm + c;
965 // select problematic cells only
966 if (hmap->GetBinContent(hmap->FindBin(absId,ri)) == 0) continue;
968 Int_t nModule, eta, phi;
969 AbsId2EtaPhi(c, nModule, eta, phi);
973 // calculate number of bad blocks
975 for (Int_t b = 0; b < 288; b++)
976 if (blk2x2[b] == 4) nbad2x2++;
979 if (nbad2x2 == 0) continue;
981 printf(" SM%i:", sm);
984 if (nbad2x2 == 288) {
985 printf(" missing the whole supermodule!\n");
990 if (nbad2x2 >= 144) {
995 for (Int_t b = 0; b < 288; b++)
996 if (blk2x2[b] == 4) nRCU[GetEMCALRCUNumber(b)]++;
998 if (nRCU[0] == 144) {
1000 for (Int_t b = 0; b < 288; b++)
1001 if (blk2x2[b] == 4 && GetEMCALRCUNumber(b) == 0) blk2x2[b] = 0;
1004 if (nRCU[1] == 144) {
1006 for (Int_t b = 0; b < 288; b++)
1007 if (blk2x2[b] == 4 && GetEMCALRCUNumber(b) == 1) blk2x2[b] = 0;
1015 for (Int_t b = 0; b < 288; b++)
1016 if (blk2x2[b] == 4) printf(" %i", b);
1017 printf(" (%i)", nbad2x2);
1022 } // supermodule loop
1028 //_________________________________________________________________________
1029 Int_t GetEMCALRCUNumber(Int_t nModule)
1031 // Returns RCU number for a 2x2 block in EMCAL;
1032 // nModule -- block number (0-287).
1034 static TH1* hRCUs = NULL;
1036 // one-time initialization
1038 hRCUs = new TH1F("hRCU1", "", 288,0,288);
1040 Int_t RCU1[144] = {8,9,10,11,20,21,22,23,32,33,34,35,44,45,46,47,56,57,58,59,68,69,70,71,80,81,82,83,
1041 92,93,94,95,104,105,106,107,116,117,118,119,128,129,130,131,140,141,142,143,148,149,150,151,
1042 152,153,154,155,160,161,162,163,164,165,166,167,172,173,174,175,176,177,178,179,184,185,186,
1043 187,188,189,190,191,196,197,198,199,200,201,202,203,208,209,210,211,212,213,214,215,220,221,
1044 222,223,224,225,226,227,232,233,234,235,236,237,238,239,244,245,246,247,248,249,250,251,256,
1045 257,258,259,260,261,262,263,268,269,270,271,272,273,274,275,280,281,282,283,284,285,286,287};
1047 for (Int_t i = 0; i < 144; i++) hRCUs->SetBinContent(RCU1[i]+1, 1);
1050 return hRCUs->GetBinContent(nModule+1);
1053 //_________________________________________________________________________
1054 void AbsId2EtaPhi(Int_t absId, Int_t &nModule, Int_t &eta, Int_t &phi, Int_t det = 0)
1056 // Converts cell absId --> (sm,eta,phi);
1058 // nModule -- 2x2 block number for EMCAL;
1059 // det -- detector: 0/EMCAL, 1/PHOS.
1063 AliEMCALGeometry *geomEMCAL = AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
1065 Int_t nSupMod, nIphi, nIeta;
1066 geomEMCAL->GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
1067 geomEMCAL->GetCellPhiEtaIndexInSModule(nSupMod, nModule, nIphi, nIeta, phi, eta);
1072 else if (det == 1) {
1073 AliPHOSGeometry *geomPHOS = AliPHOSGeometry::GetInstance("IHEP");
1076 geomPHOS->AbsToRelNumbering(absId, relid);
1087 Error("AbsId2EtaPhi", "Wrong detector");
1091 //_________________________________________________________________________
1092 void TestCellsTH1(Int_t nruns, Int_t runNumbers[], char *hname,
1093 char* htitle = "", char* hytitle = "",
1094 Int_t dnbins = 200, Double_t dmaxval = -1, Double_t goodmin = -1, Double_t goodmax = -1)
1096 // Test for bad cells by plotting a distribution of a TH1 histogram.
1097 // The histogram is obtained as a sum over runs of TH1 per run.
1099 // hname -- histogram name to process;
1100 // htitle, hytitle -- histogram title and Y axis title;
1101 // dnbins -- number of bins in distribution;
1102 // dmaxval -- X axis maximum in distribution.
1103 // goodmin,goodmax -- the region outside which a cell is considered as bad.
1105 // -1 value for maxval/goodmin/goodmax -- process a variable automatically.
1107 // initialize histogram
1108 TH1* histo = (TH1*) gInputFile->Get(Form("run%i_%s", runNumbers[0], hname));
1109 histo->SetName(hname);
1110 histo->SetTitle(htitle);
1111 histo->SetXTitle("AbsId");
1112 histo->SetYTitle(hytitle);
1115 for (Int_t i = 1; i < nruns; i++) {
1116 TH1* h = (TH1*) gInputFile->Get(Form("run%i_%s", runNumbers[i], hname));
1121 histo->Scale(1./(Double_t)GetTotalNumberOfEvents(nruns, runNumbers));
1122 Process1(histo, Form("TestCellsTH1_%s",hname), dnbins, dmaxval, goodmin, goodmax);
1125 //_________________________________________________________________________
1126 void TestCellEandN(Double_t Emin = 0.1, char* hname = "hCellAmplitude", Int_t dnbins = 200,
1127 Double_t maxval1 = -1, Double_t goodmin1 = -1, Double_t goodmax1 = -1,
1128 Double_t maxval2 = -1, Double_t goodmin2 = -1, Double_t goodmax2 = -1,
1129 Double_t maxval3 = -1, Double_t goodmin3 = -1, Double_t goodmax3 = -1)
1131 // Three more tests for bad cells:
1132 // 1) total deposited energy;
1133 // 2) total number of entries;
1134 // 3) average energy = [total deposited energy]/[total number of entries].
1136 // Based on summary histograms. Possible choises:
1137 // hCellAmplitude, hCellAmplitudeEhigh, hCellAmplitudeNonLocMax, hCellAmplitudeEhighNonLocMax
1139 // Emin -- minimum cell amplitude to count;
1140 // hname -- name (in file) of TH2 histogram with cell amplitudes;
1141 // dnbins -- number of bins in distributions;
1142 // maxval[123] -- maximum values on distributions for the criteria 1),2),3) respectively;
1143 // goodmin[123],goodmax[123] -- the regions on distributions outside those a cell is considered as bad.
1145 // input; X axis -- absId numbers
1146 TH2 *hCellAmplitude = (TH2*) gInputFile->Get(hname);
1148 // binning parameters
1149 Int_t ncells = hCellAmplitude->GetNbinsX();
1150 Double_t amin = hCellAmplitude->GetXaxis()->GetXmin();
1151 Double_t amax = hCellAmplitude->GetXaxis()->GetXmax();
1153 TH1* hCellEtotal = new TH1F(Form("%s_hCellEtotal_E%.2f",hname,Emin),
1154 Form("Total deposited energy, E > %.2f GeV",Emin), ncells,amin,amax);
1155 hCellEtotal->SetXTitle("AbsId");
1156 hCellEtotal->SetYTitle("Energy, GeV");
1158 TH1F *hCellNtotal = new TH1F(Form("%s_hCellNtotal_E%.2f",hname,Emin),
1159 Form("Total number of entries, E > %.2f GeV",Emin), ncells,amin,amax);
1160 hCellNtotal->SetXTitle("AbsId");
1161 hCellNtotal->SetYTitle("Entries");
1163 TH1F *hCellEtoNtotal = new TH1F(Form("%s_hCellEtoNtotal_E%.2f",hname,Emin),
1164 Form("Average energy per hit, E > %.2f GeV",Emin), ncells,amin,amax);
1165 hCellEtoNtotal->SetXTitle("AbsId");
1166 hCellEtoNtotal->SetYTitle("Energy, GeV");
1169 for (Int_t c = 1; c <= ncells; c++) {
1173 for (Int_t i = 1; i <= hCellAmplitude->GetNbinsY(); i++) {
1174 Double_t E = hCellAmplitude->GetYaxis()->GetBinCenter(i);
1175 Double_t N = hCellAmplitude->GetBinContent(c, i);
1176 if (E < Emin) continue;
1181 hCellEtotal->SetBinContent(c, Esum);
1182 hCellNtotal->SetBinContent(c, Nsum);
1184 if (Nsum > 0.5) // number of entries >= 1
1185 hCellEtoNtotal->SetBinContent(c, Esum/Nsum);
1188 delete hCellAmplitude;
1190 Process1(hCellEtotal, Form("%s_CellE", hname), dnbins, maxval1, goodmin1, goodmax1);
1191 Process1(hCellNtotal, Form("%s_CellN", hname), dnbins, maxval2, goodmin2, goodmax2);
1192 Process1(hCellEtoNtotal, Form("%s_CellE/N", hname), dnbins, maxval3, goodmin3, goodmax3);
1195 //_________________________________________________________________________
1196 void TestCellShapes(Double_t fitEmin, Double_t fitEmax, char* hname = "hCellAmplitude", Int_t dnbins = 1000,
1197 Double_t maxval1 = -1, Double_t goodmin1 = -1, Double_t goodmax1 = -1,
1198 Double_t maxval2 = -1, Double_t goodmin2 = -1, Double_t goodmax2 = -1,
1199 Double_t maxval3 = -1, Double_t goodmin3 = -1, Double_t goodmax3 = -1)
1201 // Test cells shape using fit function f(x)=A*exp(-B*x)/x^2.
1202 // Produce values per cell + distributions for A, B and chi2/ndf parameters.
1204 // fitEmin, fitEmax -- fit range;
1205 // hname -- name (in file) of TH2 histogram with cell amplitudes;
1206 // dnbins -- number of bins in distributions;
1207 // maxval[123] -- maximum values on distributions for the criteria A, B and chi2/ndf respectively;
1208 // goodmin[123],goodmax[123] -- the regions on distributions outside those a cell is considered as bad.
1210 // -1 value for maxval/goodmin/goodmax -- process a variable automatically.
1212 // Note: numbers are optimized for EMCAL.
1213 // TODO: check for PHOS
1215 // input; X axis -- absId numbers
1216 TH2 *hCellAmplitude = (TH2*) gInputFile->Get(hname);
1218 // binning parameters
1219 Int_t ncells = hCellAmplitude->GetNbinsX();
1220 Double_t amin = hCellAmplitude->GetXaxis()->GetXmin();
1221 Double_t amax = hCellAmplitude->GetXaxis()->GetXmax();
1223 // initialize histograms
1224 TH1 *hFitA = new TH1F(Form("hFitA_%s",hname),"Fit A value", ncells,amin,amax);
1225 hFitA->SetXTitle("AbsId");
1226 hFitA->SetYTitle("A");
1228 TH1 *hFitB = new TH1F(Form("hFitB_%s",hname),"Fit B value", ncells,amin,amax);
1229 hFitB->SetXTitle("AbsId");
1230 hFitB->SetYTitle("B");
1232 TH1 *hFitChi2Ndf = new TH1F(Form("hFitChi2Ndf_%s",hname),"Fit #chi^{2}/ndf value", ncells,amin,amax);
1233 hFitChi2Ndf->SetXTitle("AbsId");
1234 hFitChi2Ndf->SetYTitle("#chi^{2}/ndf");
1236 // total number of events; to estimate A value
1237 TH1* hNEventsProcessedPerRun = (TH1*) gInputFile->Get("hNEventsProcessedPerRun");
1238 Double_t ntotal = hNEventsProcessedPerRun->Integral(1, hNEventsProcessedPerRun->GetNbinsX());
1241 TF1 *fit = new TF1("fit", "[0]*exp(-[1]*x)/x^2");
1242 fit->SetParLimits(0, ntotal*1e-8,ntotal*1e-4);
1243 fit->SetParLimits(1, 0.,30.);
1244 fit->SetParameter(0, ntotal*1e-6);
1245 fit->SetParameter(1, 1.5);
1247 for (Int_t i = 1; i <= ncells; i++) {
1248 TH1 *hCell = hCellAmplitude->ProjectionY("",i,i);
1249 if (hCell->GetEntries() == 0) continue;
1251 hCell->Fit(fit, "0LQEM", "", fitEmin, fitEmax);
1254 hFitA->SetBinContent(i, fit->GetParameter(0));
1255 hFitB->SetBinContent(i, fit->GetParameter(1));
1256 if (fit->GetNDF() != 0)
1257 hFitChi2Ndf->SetBinContent(i, fit->GetChisquare()/fit->GetNDF());
1260 delete hCellAmplitude;
1262 // automatic parameters, if requested
1263 if (maxval1 < 0) maxval1 = 4e-6 * ntotal;
1264 if (maxval2 < 0) maxval2 = 10.;
1265 if (maxval3 < 0) maxval3 = 15.;
1267 Process1(hFitA, Form("%s_FitA", hname), dnbins, maxval1, goodmin1, goodmax1);
1268 Process1(hFitB, Form("%s_FitB", hname), dnbins, maxval2, goodmin2, goodmax2);
1269 Process1(hFitChi2Ndf, Form("%s_FitChi2ndf", hname), dnbins, maxval3, goodmin3, goodmax3);
1272 //_________________________________________________________________________
1273 void Process1(TH1* inhisto, char* label = "", Int_t dnbins = 200,
1274 Double_t dmaxval = -1, Double_t goodmin = -1, Double_t goodmax = -1)
1276 // Actual distribution analysis for a TH1 histogram:
1277 // 1) create a distribution for the input histogram;
1279 // 3) take a decision about bad cells.
1281 // inhisto -- input histogram;
1282 // label -- text label for stdout;
1283 // dnbins -- number of bins in distribution;
1284 // goodmin,goodmax -- cells outside this region are considered as bad;
1285 // dmaxval -- maximum value on distribution histogram.
1286 // The later is required in cases where a bad cell kills the whole distribution:
1287 // limiting distribution maximum value solves the problem.
1290 dmaxval = inhisto->GetMaximum()*1.01; // 1.01 - to see the last bin
1292 TH1 *distrib = new TH1F(Form("%sDistr",inhisto->GetName()), "", dnbins, inhisto->GetMinimum(), dmaxval);
1293 distrib->SetXTitle(inhisto->GetYaxis()->GetTitle());
1294 distrib->SetYTitle("Entries");
1296 // fill distribution
1297 for (Int_t c = 1; c <= inhisto->GetNbinsX(); c++)
1298 distrib->Fill(inhisto->GetBinContent(c));
1300 // draw histogram + distribution
1301 TCanvas *c1 = new TCanvas(inhisto->GetName(),inhisto->GetName(), 800,400);
1305 gPad->SetLeftMargin(0.14);
1306 gPad->SetRightMargin(0.06);
1308 inhisto->SetTitleOffset(1.7,"Y");
1312 gPad->SetLeftMargin(0.14);
1313 gPad->SetRightMargin(0.10);
1317 // simple way to estimate the left margin for good cells region:
1318 // go from left to right and find the first bin with content 2,
1319 // then go from this bin right to left while bin content is nonzero
1321 goodmin = distrib->GetXaxis()->GetXmin();
1323 for (Int_t i = 1; i <= distrib->GetNbinsX(); i++)
1324 if (distrib->GetBinContent(i) == 2) {
1325 while (i > 1 && distrib->GetBinContent(i-1) > 0) i--;
1326 goodmin = distrib->GetBinLowEdge(i);
1331 // the same automatic algorithm as above, but reflected
1333 goodmax = distrib->GetXaxis()->GetXmax();
1335 for (Int_t i = distrib->GetNbinsX(); i >= 1; i--)
1336 if (distrib->GetBinContent(i) == 2) {
1337 while (i < distrib->GetNbinsX() && distrib->GetBinContent(i+1) > 0) i++;
1338 goodmax = distrib->GetXaxis()->GetBinUpEdge(i);
1344 TLine *lline = new TLine(goodmin, 0, goodmin, distrib->GetMaximum());
1345 lline->SetLineColor(kOrange);
1346 lline->SetLineStyle(7);
1349 TLine *rline = new TLine(goodmax, 0, goodmax, distrib->GetMaximum());
1350 rline->SetLineColor(kOrange);
1351 rline->SetLineStyle(7);
1355 TLegend *leg = new TLegend(0.60,0.82,0.9,0.88);
1356 leg->AddEntry(lline, "Good region boundary","l");
1362 printf("Bad cells by criterum \"%s\":", label);
1364 // print bad cell numbers (outside goodmin,goodmax region)
1366 for (Int_t c = 1; c <= inhisto->GetNbinsX(); c++)
1367 if (inhisto->GetBinContent(c) < goodmin || inhisto->GetBinContent(c) > goodmax) {
1368 printf(" %.0f", inhisto->GetBinLowEdge(c));
1372 printf(" (%i total)\n\n", ntot);
1375 //_________________________________________________________________________
1376 void DrawCell(Int_t absId, Double_t fitEmin = 0.3, Double_t fitEmax = 1.,
1377 Double_t Emax = -1, char* hname = "hCellAmplitude")
1379 // Draw one cell spectrum with a fit.
1381 // fitEmin, fitEmax -- fit range;
1382 // Emax -- maximum value on X axis to show (-1 = no limit);
1383 // hname -- TH2 histogram name to process, possible values:
1384 // "hCellAmplitude", "hCellAmplitudeEHigh", "hCellAmplitudeNonLocMax", "hCellAmplitudeEhighNonLocMax".
1386 // input; X axis -- absId numbers
1387 TH2* hCellAmplitude = (TH2*) gInputFile->Get(hname);
1389 Int_t bin = hCellAmplitude->GetXaxis()->FindBin(absId);
1390 TH1* hCell = hCellAmplitude->ProjectionY(Form("hCell%i_%s",absId,hname),bin,bin);
1391 hCell->SetXTitle("Energy, GeV");
1392 hCell->SetYTitle("Entries");
1393 hCell->SetTitle(Form("Cell %i, %s", absId, hname));
1394 if (Emax > 0) hCell->SetAxisRange(0, Emax);
1397 TCanvas *c1 = new TCanvas(Form("hCell%i_%s",absId,hname), Form("hCell%i_%s",absId,hname), 400,400);
1398 gPad->SetLeftMargin(0.12);
1399 gPad->SetRightMargin(0.08);
1400 gPad->SetBottomMargin(0.12);
1404 // total number of events
1405 TH1* hNEventsProcessedPerRun = (TH1*) gInputFile->Get("hNEventsProcessedPerRun");
1406 Double_t ntotal = hNEventsProcessedPerRun->Integral(1, hNEventsProcessedPerRun->GetNbinsX());
1409 TF1 *fit = new TF1("fit", "[0]*exp(-[1]*x)/x^2");
1410 fit->SetLineColor(kRed);
1411 fit->SetLineWidth(2);
1412 fit->SetParName(0, "A");
1413 fit->SetParName(1, "B");
1414 fit->SetParLimits(0, ntotal*1e-8,ntotal*1e-4);
1415 fit->SetParLimits(1, 0.,30.);
1416 fit->SetParameter(0, ntotal*1e-6);
1417 fit->SetParameter(1, 1.);
1418 hCell->Fit(fit,"LQEM", "", fitEmin, fitEmax);
1421 TLegend *leg = new TLegend(0.5,0.75,0.9,0.8);
1422 leg->AddEntry(fit, "A*exp(-B*x)/x^{2}","l");
1428 //_________________________________________________________________________
1429 void DrawCellTime(Int_t absId)
1431 // Draw one cell time spectrum
1433 // input; X axis -- absId numbers
1434 TH2* hCellTime = (TH2*) gInputFile->Get("hCellTime");
1436 Int_t bin = hCellTime->GetXaxis()->FindBin(absId);
1437 TH1* hCell = hCellTime->ProjectionY(Form("hCellTime%i",absId),bin,bin);
1438 hCell->SetXTitle("Time, s");
1439 hCell->SetYTitle("Entries");
1440 hCell->SetTitle(Form("Cell %i time", absId));
1443 TCanvas *c1 = new TCanvas(Form("hCellTime%i",absId), Form("hCellTime%i",absId), 400,400);
1444 gPad->SetLeftMargin(0.12);
1445 gPad->SetRightMargin(0.10);
1446 gPad->SetBottomMargin(0.12);
1454 /* RUN AVERAGES AND RELATED FUNCTIONS
1458 //_________________________________________________________________________
1459 void DrawClusterAveragesPerRun(Int_t nruns, Int_t runNumbers[], Int_t ncellsMin = 1,
1460 Double_t eclusMin = 0.3, Double_t eclusMax = -1,
1463 // Draws cluster averages per run vs run index and number of events.
1464 // NOTE: the averages are "smoothed" a little due to a finite bin width.
1466 // ncellsMin -- minimum number of cells in cluster (>= 1);
1467 // eclusMin,eclusMax -- minimum/maximum cluster energy cut
1468 // (eclusMax = -1 means infinity);
1469 // hmap -- dead/noisy cell map, which is used for acceptance calculation due
1470 // to switched off detector parts. Acceptance is used to correct the average
1471 // number of clusters per event.
1474 TString s(Form("_NC%i_Emin=%.2fGeV",ncellsMin,eclusMin));
1475 if (eclusMax > 0) s += TString(Form("_Emax=%.2fGeV",eclusMax)).Data();
1476 if (hmap) s += "_corr4accept";
1477 char *suff = s.Data();
1479 if (eclusMax < 0) eclusMax = 1e+5;
1481 // supermodule region
1484 while (SM1 <= SM2 && !gInputFile->Get(Form("run%i_hNCellsInClusterSM%i",runNumbers[0],SM1)) ) SM1++;
1485 while (SM2 >= SM1 && !gInputFile->Get(Form("run%i_hNCellsInClusterSM%i",runNumbers[0],SM2)) ) SM2--;
1487 // initialize histograms
1488 hAvECluster = new TH1F(Form("hAvECluster%s",suff), "Average cluster energy", nruns,0,nruns);
1489 SetRunLabel(hAvECluster,nruns,runNumbers,1);
1490 // hAvECluster->SetXTitle("Run index");
1491 hAvECluster->SetTitleOffset(0.3,"Y");
1492 hAvECluster->SetYTitle("Energy, GeV");
1493 hAvECluster->SetTickLength(0.01,"Y");
1494 hAvECluster->SetLabelSize(0.06,"XY");
1495 hAvECluster->SetTitleSize(0.06,"XY");
1497 hAvNCluster = new TH1F(Form("hAvNCluster%s",suff), "Average number of clusters per event", nruns,0,nruns);
1498 SetRunLabel(hAvNCluster,nruns,runNumbers,1);
1499 // hAvNCluster->SetXTitle("Run index");
1500 hAvNCluster->SetTitleOffset(0.3,"Y");
1501 hAvNCluster->SetYTitle("Number of clusters");
1502 hAvNCluster->SetTickLength(0.01,"Y");
1503 hAvNCluster->SetLabelSize(0.06,"XY");
1504 hAvNCluster->SetTitleSize(0.06,"XY");
1506 hAvNCellsInCluster = new TH1F(Form("hAvNCellsInCluster%s",suff), "Average number of cells in cluster", nruns,0,nruns);
1507 SetRunLabel(hAvNCellsInCluster,nruns,runNumbers,1);
1508 // hAvNCellsInCluster->SetXTitle("Run index");
1509 hAvNCellsInCluster->SetTitleOffset(0.3,"Y");
1510 hAvNCellsInCluster->SetYTitle("Number of cells");
1511 hAvNCellsInCluster->SetTickLength(0.01,"Y");
1512 hAvNCellsInCluster->SetLabelSize(0.06,"XY");
1513 hAvNCellsInCluster->SetTitleSize(0.06,"XY");
1515 // initialize per SM histograms
1516 TH1* hAvEClusterSM[10];
1517 TH1* hAvNClusterSM[10];
1518 TH1* hAvNCellsInClusterSM[10];
1520 for (Int_t sm = SM1; sm <= SM2; sm++) {
1521 hAvEClusterSM[sm] = new TH1F(Form("hAvEClusterSM%i%s",sm,suff),"", nruns,0,nruns);
1522 hAvNClusterSM[sm] = new TH1F(Form("hAvNClusterSM%i%s",sm,suff),"", nruns,0,nruns);
1523 hAvNCellsInClusterSM[sm] = new TH1F(Form("hAvNCellsInClusterSM%i%s",sm,suff),"", nruns,0,nruns);
1526 // fill all the histograms per run index
1527 for (Int_t ri = 0; ri < nruns; ri++)
1529 Int_t nevents = GetNumberOfEvents(runNumbers[ri]);
1531 // number of switched off supermodules
1534 Double_t Eclus_total = 0; // total cluster energy
1535 Double_t Nclus_total = 0; // total number of clusters
1536 Double_t Ncells_total = 0; // total number of cells
1539 for (Int_t sm = SM1; sm <= SM2; sm++) {
1540 TH2* hNCellsInClusterSM = (TH2*) gInputFile->Get(Form("run%i_hNCellsInClusterSM%i",runNumbers[ri],sm));
1541 if (hNCellsInClusterSM == 0) {
1542 Error("DrawClusterAveragesPerRun","Run %d does contain histogram %s",runNumbers[ri],Form("run%i_hNCellsInClusterSM%i",runNumbers[ri],sm));
1546 // the same as above, but per supermodule
1547 Double_t Eclus_totalSM = 0;
1548 Double_t Nclus_totalSM = 0;
1549 Double_t Ncells_totalSM = 0;
1551 // X axis -- cluster energy, Y axis -- number of cells
1552 for (Int_t x = 1; x <= hNCellsInClusterSM->GetNbinsX(); x++)
1553 for (Int_t y = 1+ncellsMin; y <= hNCellsInClusterSM->GetNbinsY(); y++) {//NOTE: bin 1 correspond to ncellsMin=0
1554 Double_t Eclus = hNCellsInClusterSM->GetXaxis()->GetBinCenter(x);
1555 Double_t Ncells = hNCellsInClusterSM->GetYaxis()->GetBinLowEdge(y);
1556 Double_t Nclus = hNCellsInClusterSM->GetBinContent(x,y);
1558 // cut on cluster energy
1559 if (Eclus < eclusMin || Eclus > eclusMax) continue;
1561 Eclus_totalSM += Eclus * Nclus;
1562 Nclus_totalSM += Nclus;
1563 Ncells_totalSM += Ncells * Nclus;
1566 delete hNCellsInClusterSM;
1568 // correct for acceptance
1570 Double_t accep = GetAcceptance(sm, hmap, ri);
1572 Eclus_totalSM /= accep;
1573 Nclus_totalSM /= accep;
1574 Ncells_totalSM /= accep;
1579 Eclus_total += Eclus_totalSM;
1580 Nclus_total += Nclus_totalSM;
1581 Ncells_total += Ncells_totalSM;
1583 hAvNClusterSM[sm]->SetBinContent(ri+1, Nclus_totalSM/nevents);
1584 if (Nclus_totalSM > 0) hAvEClusterSM[sm]->SetBinContent(ri+1, Eclus_totalSM/Nclus_totalSM);
1585 if (Nclus_totalSM > 0) hAvNCellsInClusterSM[sm]->SetBinContent(ri+1, Ncells_totalSM/Nclus_totalSM);
1586 } // supermodule loop
1588 hAvNCluster->SetBinContent(ri+1, Nclus_total/nevents/(SM2-SM1+1-noSM));
1589 if (Nclus_total > 0) hAvECluster->SetBinContent(ri+1, Eclus_total/Nclus_total);
1590 if (Nclus_total > 0) hAvNCellsInCluster->SetBinContent(ri+1, Ncells_total/Nclus_total);
1594 /* Draw results vs run index
1597 TCanvas *c1 = new TCanvas(Form("ClusterAveragesRuns%s",suff),
1598 Form("ClusterAveragesRuns%s",suff), 1000,707);
1601 // average cluster energy
1603 gPad->SetLeftMargin(0.04);
1604 gPad->SetRightMargin(0.02);
1605 gPad->SetTopMargin(0.10);
1606 gPad->SetBottomMargin(0.14);
1609 TLegend *leg = new TLegend(0.65,0.16,0.75,0.15+0.08*(SM2-SM1+1));
1611 hAvECluster->SetAxisRange(hAvECluster->GetMinimum()*0.5, hAvECluster->GetMaximum()*1.5,"Y");
1612 hAvECluster->SetLineWidth(2);
1613 hAvECluster->Draw();
1614 leg->AddEntry(hAvECluster, "All Modules","l");
1615 for (Int_t sm = SM1; sm <= SM2; sm++) {
1616 hAvEClusterSM[sm]->SetLineColor(colorsSM[sm]);
1617 hAvEClusterSM[sm]->SetLineWidth(1);
1618 hAvEClusterSM[sm]->Draw("same");
1619 leg->AddEntry(hAvEClusterSM[sm], Form("Module %i",5-sm),"l");
1621 hAvECluster->Draw("same"); // to top
1624 // average number of clusters
1626 gPad->SetLeftMargin(0.04);
1627 gPad->SetRightMargin(0.02);
1628 gPad->SetTopMargin(0.10);
1629 gPad->SetBottomMargin(0.14);
1632 leg = new TLegend(0.65,0.16,0.75,0.15+0.08*(SM2-SM1+1));
1634 hAvNCluster->SetAxisRange(0, hAvNCluster->GetMaximum()*1.5,"Y");
1635 hAvNCluster->SetLineWidth(2);
1636 hAvNCluster->Draw();
1637 leg->AddEntry(hAvNCluster, Form("(All Moduls)/%i",SM2-SM1+1),"l");
1638 for (Int_t sm = SM1; sm <= SM2; sm++) {
1639 hAvNClusterSM[sm]->SetLineColor(colorsSM[sm]);
1640 hAvNClusterSM[sm]->SetLineWidth(1);
1641 hAvNClusterSM[sm]->Draw("same");
1642 leg->AddEntry(hAvNClusterSM[sm], Form("Module %i",5-sm),"l");
1644 hAvNCluster->Draw("same"); // to top
1647 // average number of cells in cluster
1649 gPad->SetLeftMargin(0.04);
1650 gPad->SetRightMargin(0.02);
1651 gPad->SetTopMargin(0.10);
1652 gPad->SetBottomMargin(0.14);
1655 leg = new TLegend(0.65,0.16,0.75,0.15+0.08*(SM2-SM1+1));
1657 hAvNCellsInCluster->SetAxisRange(0, hAvNCellsInCluster->GetMaximum()*1.3,"Y");
1658 hAvNCellsInCluster->SetLineWidth(2);
1659 hAvNCellsInCluster->Draw();
1660 leg->AddEntry(hAvNCellsInCluster, "All Modules","l");
1661 for (Int_t sm = SM1; sm <= SM2; sm++) {
1662 hAvNCellsInClusterSM[sm]->SetLineColor(colorsSM[sm]);
1663 hAvNCellsInClusterSM[sm]->SetLineWidth(1);
1664 hAvNCellsInClusterSM[sm]->Draw("same");
1665 leg->AddEntry(hAvNCellsInClusterSM[sm], Form("Module %i",5-sm),"l");
1667 hAvNCellsInCluster->Draw("same"); // to top
1671 if (!trendPlotEvents) return;
1673 /* Draw the same vs number of events
1676 TCanvas *c2 = new TCanvas(Form("ClusterAveragesEvents%s",suff),
1677 Form("ClusterAveragesEvents%s",suff), 1000,707);
1680 // average cluster energy
1682 gPad->SetLeftMargin(0.04);
1683 gPad->SetRightMargin(0.02);
1684 gPad->SetTopMargin(0.10);
1685 gPad->SetBottomMargin(0.13);
1686 leg = new TLegend(0.65,0.15,0.85,0.15+0.04*(SM2-SM1+1));
1688 TH1* hAvEClusterEvents = ConvertMapRuns2Events(nruns, runNumbers, hAvECluster);
1689 hAvEClusterEvents->Draw();
1690 leg->AddEntry(hAvEClusterEvents, "All Modules","l");
1691 for (Int_t sm = SM1; sm <= SM2; sm++) {
1692 TH1* hAvEClusterSMEvents = ConvertMapRuns2Events(nruns, runNumbers, hAvEClusterSM[sm]);
1693 hAvEClusterSMEvents->Draw("same");
1694 leg->AddEntry(hAvEClusterSMEvents, Form("Module %i",5-sm),"l");
1696 hAvEClusterEvents->Draw("same"); // to top
1699 // average number of clusters
1701 gPad->SetLeftMargin(0.04);
1702 gPad->SetRightMargin(0.02);
1703 gPad->SetTopMargin(0.10);
1704 gPad->SetBottomMargin(0.13);
1705 leg = new TLegend(0.65,0.15,0.85,0.15+0.04*(SM2-SM1+1));
1707 TH1* hAvNClusterEvents = ConvertMapRuns2Events(nruns, runNumbers, hAvNCluster);
1708 hAvNClusterEvents->Draw();
1709 leg->AddEntry(hAvNClusterEvents, Form("(All Modules)/%i",SM2-SM1+1),"l");
1710 for (Int_t sm = SM1; sm <= SM2; sm++) {
1711 TH1* hAvNClusterSMEvents = ConvertMapRuns2Events(nruns, runNumbers, hAvNClusterSM[sm]);
1712 hAvNClusterSMEvents->Draw("same");
1713 leg->AddEntry(hAvNClusterSMEvents, Form("Module %i",5-sm),"l");
1715 hAvNClusterEvents->Draw("same"); // to top
1718 // average number of cells in cluster
1720 gPad->SetLeftMargin(0.04);
1721 gPad->SetRightMargin(0.02);
1722 gPad->SetTopMargin(0.10);
1723 gPad->SetBottomMargin(0.13);
1724 leg = new TLegend(0.65,0.15,0.85,0.15+0.04*(SM2-SM1+1));
1726 TH1* hAvNCellsInClusterEvents = ConvertMapRuns2Events(nruns, runNumbers, hAvNCellsInCluster);
1727 hAvNCellsInClusterEvents->Draw();
1728 leg->AddEntry(hAvNCellsInClusterEvents, "All Modules","l");
1729 for (Int_t sm = SM1; sm <= SM2; sm++) {
1730 TH1* hAvNCellsInClusterSMEvents = ConvertMapRuns2Events(nruns, runNumbers, hAvNCellsInClusterSM[sm]);
1731 hAvNCellsInClusterSMEvents->Draw("same");
1732 leg->AddEntry(hAvNCellsInClusterSMEvents, Form("Module %i",5-sm),"l");
1734 hAvNCellsInClusterEvents->Draw("same"); // to top
1740 //_________________________________________________________________________
1741 Double_t GetAcceptance(Int_t sm, TH2* hmap, Int_t ri)
1743 // Returns [#cells - #dead]/#cells for a supermodule.
1744 // hmap -- dead/noisy cell map;
1747 // guess number of cells per supermodule
1748 Int_t nSM = 1152; // EMCAL
1749 if (hmap->GetXaxis()->GetXmin() == 1) {// PHOS
1751 sm--; // starts from 1, convenient from 0
1756 for (Int_t k = 1; k <= nSM; k++)
1757 if (hmap->GetBinContent(nSM*sm + k, ri+1) < 0)
1760 return ((Double_t) (nSM - ndead))/nSM;
1763 //_________________________________________________________________________
1764 void DrawPi0Slice(Int_t run, Int_t sm = -1)
1766 // Draw the pi0 peak;
1767 // run,sm -- run number and supermodule to take;
1768 // sm < 0 -- draw for whole detector.
1771 if (sm >= 0) {//particular supermodule
1772 h = (TH1*) gInputFile->Get(Form("run%i_hPi0MassSM%iSM%i",run,sm,sm));
1774 Error("DrawPi0Slice","Run %d does contain histogram %s",run,Form("run%i_hPi0MassSM%iSM%i",run,sm,sm));
1777 h->SetName(Form("hPi0SliceSM%i_run%i",sm,run));
1778 h->SetTitle(Form("#pi^{0} in M%i, run %i, %.0fk events", 5-sm, run, GetNumberOfEvents(run)/1e+3));
1780 else {// whole detector
1781 for (Int_t sm1 = 0; sm1 < 10; sm1++)
1782 for (Int_t sm2 = sm1; sm2 < 10; sm2++) {
1783 TH1* one = (TH1*) gInputFile->Get(Form("run%i_hPi0MassSM%iSM%i",run,sm1,sm2));
1792 h->SetName(Form("hPi0Slice_run%i",run));
1793 h->SetTitle(Form("#pi^{0} in all modules, run %i, %.0fk events", run, GetNumberOfEvents(run)/1e+3));
1796 h->SetXTitle("M_{#gamma#gamma}, GeV");
1797 h->SetYTitle("Entries");
1798 h->SetTitleOffset(1.7,"Y");
1801 if (sm >= 0) c1 = new TCanvas(Form("hPi0SliceSM%i_run%i",sm,run),Form("hPi0SliceSM%i_run%i",sm,run), 400,400);
1802 else c1 = new TCanvas(Form("hPi0Slice_run%i",run),Form("hPi0Slice_run%i",run), 400,400);
1804 gPad->SetLeftMargin(0.14);
1805 gPad->SetRightMargin(0.06);
1808 Double_t nraw, enraw, mass, emass, sigma, esigma;
1809 FitPi0(h, nraw, enraw, mass, emass, sigma, esigma);
1812 TF1* fitfun = h->GetFunction("fitfun");
1814 Double_t emin, emax;
1815 fitfun->GetRange(emin, emax);
1817 backgr = new TF1("mypol2", "[0] + [1]*(x-0.135) + [2]*(x-0.135)^2", emin, emax);
1818 backgr->SetLineColor(kBlue);
1819 backgr->SetLineWidth(2);
1820 backgr->SetLineStyle(3);
1821 backgr->SetParameters(fitfun->GetParameter(3), fitfun->GetParameter(4), fitfun->GetParameter(5));
1822 backgr->Draw("same");
1828 //_________________________________________________________________________
1829 void DrawPi0Averages(Int_t nruns, Int_t runNumbers[], Bool_t samesm = kFALSE, TH2* hmap = NULL)
1831 // Draw average number of pi0s per event, pi0 mass position and pi0 width per run index.
1832 // Errors are also drawn.
1834 // samesm -- take only pi0s for which gammas were is same SM;
1835 // hmap -- if not NULL, do simple (area law based) acceptance correction.
1837 // TODO: PHOS needs pi0 between SMs rather than in one SM
1840 char *suff = TString(Form("%s%s", samesm ? "_sameSM":"", hmap ? "_corr4accep":"")).Data();
1842 // supermodule region
1845 while (SM1 <= SM2 && !gInputFile->Get(Form("run%i_hPi0MassSM%iSM%i",runNumbers[0],SM1,SM1)) ) SM1++;
1846 while (SM2 >= SM1 && !gInputFile->Get(Form("run%i_hPi0MassSM%iSM%i",runNumbers[0],SM2,SM2)) ) SM2--;
1848 // initialize histograms for the entire detector
1849 TH1* hPi0Num = new TH1F(Form("hPi0Num%s",suff),"Average number of #pi^{0}s per event", nruns,0,nruns);
1850 SetRunLabel(hPi0Num,nruns,runNumbers,1);
1851 // hPi0Num->SetXTitle("Run index");
1852 hPi0Num->SetYTitle("Number of #pi^{0}s");
1854 TH1* hPi0Mass = new TH1F(Form("hPi0Mass%s",suff),"#pi^{0} mass position", nruns,0,nruns);
1855 SetRunLabel(hPi0Mass,nruns,runNumbers,1);
1856 // hPi0Mass->SetXTitle("Run index");
1857 hPi0Mass->SetYTitle("M_{#pi^{0}}, GeV/c^{2}");
1859 TH1* hPi0Sigma = new TH1F(Form("hPi0Sigma%s",suff),"#pi^{0} width", nruns,0,nruns);
1860 SetRunLabel(hPi0Sigma,nruns,runNumbers,1);
1861 // hPi0Sigma->SetXTitle("Run index");
1862 hPi0Sigma->SetYTitle("#sigma_{#pi^{0}}, GeV/c^{2}");
1864 // initialize histograms per SM
1866 TH1* hPi0MassSM[10];
1867 TH1* hPi0SigmaSM[10];
1869 for (Int_t sm = SM1; sm <= SM2; sm++) {
1870 hPi0NumSM[sm] = new TH1F(Form("hPi0NumSM%i%s",sm,suff),"", nruns,0,nruns);
1871 hPi0MassSM[sm] = new TH1F(Form("hPi0MassSM%i%s",sm,suff),"", nruns,0,nruns);
1872 hPi0SigmaSM[sm] = new TH1F(Form("hPi0SigmaSM%i%s",sm,suff),"", nruns,0,nruns);
1875 TCanvas cFit("cFit","cFit");
1877 for (Int_t ri = 0; ri < nruns; ri++)
1879 fprintf(stderr, "\rDrawPi0Averages(): analysing run index %i/%i, run # %i... ", ri, nruns-1,runNumbers[ri]);
1881 Int_t nevents = GetNumberOfEvents(runNumbers[ri]);
1884 for (Int_t sm = SM1; sm <= SM2; sm++) {
1885 TH1* h = (TH1*) gInputFile->Get(Form("run%i_hPi0MassSM%iSM%i",runNumbers[ri],sm,sm));
1887 Error("DrawPi0Averages","Run %d does contain histogram %s",runNumbers[ri],Form("run%i_hPi0MassSM%iSM%i",runNumbers[ri],sm,sm));
1891 // supermodule acceptance
1892 Double_t accep = 1.;
1893 if (hmap) accep = GetAcceptance(sm, hmap, ri);
1894 if (accep == 0) continue; // missing SM
1896 Double_t nraw, enraw, mass, emass, sigma, esigma;
1897 FitPi0(h, nraw, enraw, mass, emass, sigma, esigma);
1899 hPi0NumSM[sm] ->SetBinContent(ri+1, nraw/nevents/accep);
1900 hPi0MassSM[sm] ->SetBinContent(ri+1, mass);
1901 hPi0SigmaSM[sm]->SetBinContent(ri+1, sigma);
1903 hPi0NumSM[sm] ->SetBinError(ri+1, enraw/nevents/accep);
1904 hPi0MassSM[sm] ->SetBinError(ri+1, emass);
1905 hPi0SigmaSM[sm]->SetBinError(ri+1, esigma);
1908 } // supermodule loop
1911 /* fill for the entire detector
1913 TH1* hsum = (TH1*) gInputFile->Get(Form("run%i_hPi0MassSM%iSM%i",runNumbers[ri],SM1,SM1));
1915 Error("DrawPi0Averages","Run %d does contain histogram %s",runNumbers[ri],Form("run%i_hPi0MassSM%iSM%i",runNumbers[ri],SM1,SM1));
1918 hsum->SetName("hSumTMP");
1920 // for acceptance correction
1923 for (Int_t sm1 = SM1; sm1 <= SM2; sm1++)
1924 for (Int_t sm2 = sm1; sm2 <= SM2; sm2++) {
1925 if (sm1 == SM1 && sm2 == SM2) continue;
1926 if (samesm && sm1 != sm2) continue;
1928 TH1* h = (TH1*) gInputFile->Get(Form("run%i_hPi0MassSM%iSM%i",runNumbers[ri],sm1,sm2));
1930 Error("DrawPi0Averages","Run %d does contain histogram %s",runNumbers[ri],Form("run%i_hPi0MassSM%iSM%i",runNumbers[ri],SM1,SM1));
1934 // correct for SM acceptance
1936 Double_t accep = GetAcceptance(sm1, hmap, ri);
1937 if (accep > 0) h->Scale(1./accep);
1945 Double_t nraw, enraw, mass, emass, sigma, esigma;
1946 FitPi0(hsum, nraw, enraw, mass, emass, sigma, esigma);
1948 hPi0Num ->SetBinContent(ri+1, nraw/nevents/(SM2-SM1+1-noSM));
1949 hPi0Mass ->SetBinContent(ri+1, mass);
1950 hPi0Sigma->SetBinContent(ri+1, sigma);
1952 hPi0Num ->SetBinError(ri+1, enraw/nevents/(SM2-SM1+1-noSM));
1953 hPi0Mass ->SetBinError(ri+1, emass);
1954 hPi0Sigma->SetBinError(ri+1, esigma);
1959 fprintf(stderr, "\n");
1965 // number of pi0s vs run index
1966 TCanvas *c1 = new TCanvas(Form("hPi0NumRuns%s",suff),Form("hPi0NumRuns%s",suff), 800,400);
1967 gPad->SetLeftMargin(0.08);
1968 gPad->SetRightMargin(0.03);
1969 gPad->SetTopMargin(0.12);
1970 gPad->SetBottomMargin(0.12);
1973 // TLegend *leg = new TLegend(0.75,0.15,0.89,0.15+0.06*(SM2-SM1+1));
1974 TLegend *leg = new TLegend(0.75,0.89-0.06*(SM2-SM1+1),0.89,0.89);
1976 hPi0Num->SetAxisRange(0., hPi0Num->GetMaximum()*1.8, "Y");
1977 hPi0Num->SetTitleOffset(0.8, "Y");
1978 hPi0Num->SetLineWidth(2);
1979 hPi0Num->SetLabelSize(0.04,"X");
1981 leg->AddEntry(hPi0Num, Form("(All Modules)/%i",SM2-SM1+1),"l");
1982 for (Int_t sm = SM1; sm <= SM2; sm++) {
1983 hPi0NumSM[sm]->SetLineColor(colorsSM[sm]);
1984 hPi0NumSM[sm]->SetLineWidth(2);
1985 hPi0NumSM[sm]->Draw("same");
1986 leg->AddEntry(hPi0NumSM[sm], Form("Module %i",5-sm),"l");
1988 hPi0Num->Draw("same"); // to the top
1989 hPi0Num->Draw("hist,same");
1993 if (trendPlotEvents) {
1995 // number of pi0s vs event count
1996 TCanvas *c2 = new TCanvas(Form("hPi0NumEvents%s",suff),Form("hPi0NumEvents%s",suff), 800,400);
1997 gPad->SetLeftMargin(0.08);
1998 gPad->SetRightMargin(0.06);
1999 gPad->SetTopMargin(0.12);
2000 leg = new TLegend(0.75,0.15,0.92,0.15+0.04*(SM2-SM1+1));
2002 TH1* hPi0NumEvents = ConvertMapRuns2Events(nruns, runNumbers, hPi0Num);
2003 hPi0NumEvents->Draw();
2004 leg->AddEntry(hPi0NumEvents, Form("(All Modules)/%i",SM2-SM1+1),"l");
2005 for (Int_t sm = SM1; sm <= SM2; sm++) {
2006 TH1* hPi0NumSMEvents = ConvertMapRuns2Events(nruns, runNumbers, hPi0NumSM[sm]);
2007 hPi0NumSMEvents->Draw("same");
2008 leg->AddEntry(hPi0NumSMEvents, Form("Module %i",5-sm),"l");
2010 hPi0NumEvents->Draw("same"); // to the top
2011 hPi0NumEvents->Draw("hist,same");
2016 // pi0 mass and width vs run index
2017 TCanvas *c3 = new TCanvas(Form("hPi0MassSigmaRuns%s",suff),Form("hPi0MassSigmaRuns%s",suff), 1000,707);
2021 gPad->SetLeftMargin(0.06);
2022 gPad->SetRightMargin(0.04);
2023 gPad->SetTopMargin(0.10);
2024 gPad->SetBottomMargin(0.13);
2027 // leg = new TLegend(0.75,0.15,0.85,0.16+0.06*(SM2-SM1+1));
2028 leg = new TLegend(0.75,0.89-0.06*(SM2-SM1+1),0.85,0.89);
2030 hPi0Mass->SetAxisRange(0.125, 0.149, "Y");
2031 hPi0Mass->SetTitleOffset(0.8, "Y");
2032 hPi0Mass->SetLineWidth(2);
2033 hPi0Mass->SetLabelSize(0.06,"X");
2035 leg->AddEntry(hPi0Mass, "All Modules","l");
2036 for (Int_t sm = SM1; sm <= SM2; sm++) {
2037 hPi0MassSM[sm]->SetLineColor(colorsSM[sm]);
2038 hPi0MassSM[sm]->SetLineWidth(2);
2039 hPi0MassSM[sm]->Draw("same");
2040 leg->AddEntry(hPi0MassSM[sm], Form("Module %i",5-sm),"l");
2042 hPi0Mass->Draw("same"); // to the top
2043 hPi0Mass->Draw("hist,same");
2047 gPad->SetLeftMargin(0.06);
2048 gPad->SetRightMargin(0.04);
2049 gPad->SetTopMargin(0.10);
2050 gPad->SetBottomMargin(0.13);
2053 // leg = new TLegend(0.75,0.15,0.85,0.16+0.06*(SM2-SM1+1));
2054 leg = new TLegend(0.75,0.89-0.06*(SM2-SM1+1),0.85,0.89);
2056 hPi0Sigma->SetAxisRange(0., hPi0Sigma->GetMaximum()*1.5, "Y");
2057 hPi0Sigma->SetTitleOffset(0.8, "Y");
2058 hPi0Sigma->SetLineWidth(2);
2059 hPi0Sigma->SetLabelSize(0.06,"X");
2061 leg->AddEntry(hPi0Sigma, "All Modules","l");
2062 for (Int_t sm = SM1; sm <= SM2; sm++) {
2063 hPi0SigmaSM[sm]->SetLineColor(colorsSM[sm]);
2064 hPi0SigmaSM[sm]->SetLineWidth(2);
2065 hPi0SigmaSM[sm]->Draw("same");
2066 leg->AddEntry(hPi0SigmaSM[sm], Form("Module %i",5-sm),"l");
2068 hPi0Sigma->Draw("same"); // to the top
2069 hPi0Sigma->Draw("hist,same");
2074 if (trendPlotEvents) {
2075 // pi0 mass and width vs number of events
2076 TCanvas *c4 = new TCanvas(Form("hPi0MassSigmaEvents%s",suff),Form("hPi0MassSigmaEvents%s",suff), 800,400);
2080 gPad->SetLeftMargin(0.16);
2081 gPad->SetRightMargin(0.08);
2082 leg = new TLegend(0.75,0.15,0.91,0.15+0.04*(SM2-SM1+1));
2084 TH1* hPi0MassEvents = ConvertMapRuns2Events(nruns, runNumbers, hPi0Mass);
2085 hPi0MassEvents->Draw();
2086 leg->AddEntry(hPi0MassEvents, "All Modules","l");
2087 for (Int_t sm = SM1; sm <= SM2; sm++) {
2088 TH1* hPi0MassSMEvents = ConvertMapRuns2Events(nruns, runNumbers, hPi0MassSM[sm]);
2089 hPi0MassSMEvents->Draw("same");
2090 leg->AddEntry(hPi0MassSMEvents, Form("Module %i",5-sm),"l");
2092 hPi0MassEvents->Draw("same"); // to the top
2093 hPi0MassEvents->Draw("hist,same");
2097 gPad->SetLeftMargin(0.16);
2098 gPad->SetRightMargin(0.08);
2099 leg = new TLegend(0.75,0.15,0.91,0.15+0.04*(SM2-SM1+1));
2101 TH1* hPi0SigmaEvents = ConvertMapRuns2Events(nruns, runNumbers, hPi0Sigma);
2102 hPi0SigmaEvents->Draw();
2103 leg->AddEntry(hPi0SigmaEvents, "All Modules","l");
2104 for (Int_t sm = SM1; sm <= SM2; sm++) {
2105 TH1* hPi0SigmaSMEvents = ConvertMapRuns2Events(nruns, runNumbers, hPi0SigmaSM[sm]);
2106 hPi0SigmaSMEvents->Draw("same");
2107 leg->AddEntry(hPi0SigmaSMEvents, Form("Module %i",5-sm),"l");
2109 hPi0SigmaEvents->Draw("same"); // to the top
2110 hPi0SigmaEvents->Draw("hist,same");
2117 //_________________________________________________________________________
2118 void DrawPi0Distributions(char *suff, Int_t nbins = 100)
2120 // Draw distributions for
2121 // 1) average number of pi0s per event;
2122 // 2) pi0 mass position;
2125 // Must be called after DrawPi0Averages() because it
2126 // searches for the corresponding histograms by name.
2128 // suff -- histograms suffix;
2129 // nbins -- number of bins in distributions.
2131 TCanvas *c1 = new TCanvas(Form("Pi0Distributions%s",suff),Form("Pi0Distributions%s",suff), 999,333);
2135 c1->cd(1)->SetLogy();
2136 gPad->SetLeftMargin(0.16);
2137 gPad->SetRightMargin(0.04);
2138 TH1* hPi0Num = (TH1*) gROOT->FindObject(Form("hPi0Num%s",suff));
2139 MakeDistribution(hPi0Num,nbins)->Draw();
2142 c1->cd(2)->SetLogy();
2143 gPad->SetLeftMargin(0.16);
2144 gPad->SetRightMargin(0.04);
2145 TH1* hPi0Mass = (TH1*) gROOT->FindObject(Form("hPi0Mass%s",suff));
2146 MakeDistribution(hPi0Mass,nbins)->Draw();
2149 c1->cd(3)->SetLogy();
2150 gPad->SetLeftMargin(0.16);
2151 gPad->SetRightMargin(0.04);
2152 TH1* hPi0Sigma = (TH1*) gROOT->FindObject(Form("hPi0Sigma%s",suff));
2153 MakeDistribution(hPi0Sigma,nbins)->Draw();
2158 //_________________________________________________________________________
2159 TH1* MakeDistribution(TH1* inhisto, Int_t nbins)
2161 // Returns distribution for the input histogram.
2163 // inhisto -- input histogram;
2164 // nbins -- number of bins in distribution.
2166 // initialize distribution; 1.01 -- to see the last bin
2167 TH1* distr = new TH1F(Form("%sDistr",inhisto->GetName()),"",
2168 nbins, inhisto->GetMinimum(),inhisto->GetMaximum()*1.01);
2169 distr->SetXTitle(inhisto->GetYaxis()->GetTitle());
2170 distr->SetYTitle("Number of runs");
2172 // fill distribution
2173 for (Int_t i = 1; i <= inhisto->GetNbinsX(); i++)
2174 distr->Fill(inhisto->GetBinContent(i));
2179 //_________________________________________________________________________
2180 void FitPi0(TH1* h, Double_t &nraw, Double_t &enraw,
2181 Double_t &mass, Double_t &emass,
2182 Double_t &sigma, Double_t &esigma,
2183 Double_t emin = 0.05, Double_t emax = 0.3, Int_t rebin = 1)
2185 // Fits the pi0 peak with crystal ball + pol2,
2186 // fills number of pi0s, mass, width and their errors.
2192 if (h->GetEntries() == 0) return NULL;
2194 if (rebin > 1) h->Rebin(rebin);
2196 // crystal ball parameters (fixed by hand for EMCAL); TODO: PHOS parameters?
2197 Double_t alpha = 1.1; // alpha >= 0
2198 Double_t n = 2.; // n > 1
2200 // CB tail parameters
2201 Double_t a = pow((n/alpha), n) * TMath::Exp(-alpha*alpha/2.);
2202 Double_t b = n/alpha - alpha;
2204 // integral value under crystal ball with amplitude = 1, sigma = 1
2205 // (will be sqrt(2pi) at alpha = infinity)
2206 Double_t nraw11 = a * pow(b+alpha, 1.-n)/(n-1.) + TMath::Sqrt(TMath::Pi()/2.) * TMath::Erfc(-alpha/TMath::Sqrt(2.));
2208 // signal (crystal ball);
2209 new TF1("cball", Form("(x-[1])/[2] > -%f ? \
2210 [0]*exp(-(x-[1])*(x-[1])/(2*[2]*[2])) \
2211 : [0]*%f*(%f-(x-[1])/[2])^(-%f)", alpha, a, b, n));
2214 new TF1("mypol2", "[0] + [1]*(x-0.135) + [2]*(x-0.135)^2", emin, emax);
2216 // signal + background
2217 TF1 *fitfun = new TF1("fitfun", "cball + mypol2", emin, emax);
2218 fitfun->SetParNames("A", "M", "#sigma", "a_{0}", "a_{1}", "a_{2}");
2219 fitfun->SetLineColor(kRed);
2220 fitfun->SetLineWidth(2);
2222 // make a preliminary fit to estimate parameters
2223 TF1* ff = new TF1("fastfit", "gaus(0) + [3]");
2224 ff->SetParLimits(0, 0., h->GetMaximum()*1.5);
2225 ff->SetParLimits(1, 0.1, 0.2);
2226 ff->SetParLimits(2, 0.004,0.030);
2227 ff->SetParameters(h->GetMaximum()/3., 0.135, 0.010, 0.);
2228 h->Fit(ff, "0QL", "", 0.105, 0.165);
2230 fitfun->SetParLimits(0, 0., h->GetMaximum()*1.5);
2231 fitfun->SetParLimits(1, 0.12, 0.15);
2232 fitfun->SetParLimits(2, 0.004,0.030);
2233 fitfun->SetParameters(ff->GetParameter(0), ff->GetParameter(1), ff->GetParameter(2), ff->GetParameter(3));
2234 h->Fit(fitfun,"QLR", "");
2236 // calculate pi0 values
2237 mass = fitfun->GetParameter(1);
2238 emass = fitfun->GetParError(1);
2240 sigma = fitfun->GetParameter(2);
2241 esigma = fitfun->GetParError(2);
2243 Double_t A = fitfun->GetParameter(0);
2244 Double_t eA = fitfun->GetParError(0);
2246 nraw = nraw11 * A * sigma / h->GetBinWidth(1);
2247 enraw = nraw * (eA/A + esigma/sigma);
2251 /* COMMON FUNCTIONS FOR RUN NUMBERS EXTRACTION, VISUALIZATION AND FILTERING.
2253 * Input: histogram hNEventsProcessedPerRun which is taken from gInputFile.
2256 //_________________________________________________________________________
2257 void GetRunNumbers(Int_t &nruns, Int_t runNumbers[])
2259 // Fill runNumbers array with run numbers according to hNEventsProcessedPerRun
2260 // histogram: actually analysed runs have positive bin content.
2262 TH1* hNEventsProcessedPerRun = (TH1*) gInputFile->Get("hNEventsProcessedPerRun");
2264 // check underflow/overflow
2265 if (hNEventsProcessedPerRun->GetBinContent(0) > 0.5)
2266 Warning("GetRunNumbers", "some run numbers are in underflow; they will be absent in the list of runs");
2267 if (hNEventsProcessedPerRun->GetBinContent(hNEventsProcessedPerRun->GetNbinsX()+1) > 0.5)
2268 Warning("GetRunNumbers", "some run numbers are in overflow; they will be absent in the list of runs");
2272 for (Int_t i = 1; i <= hNEventsProcessedPerRun->GetNbinsX(); i++)
2273 if (hNEventsProcessedPerRun->GetBinContent(i) > 0.5) {
2274 runNumbers[nruns] = hNEventsProcessedPerRun->GetBinLowEdge(i);
2279 //_________________________________________________________________________
2280 Int_t GetNumberOfEvents(Int_t run)
2282 // Return number of events in run;
2283 // run -- run number.
2285 TH1* hNEventsProcessedPerRun = (TH1*) gInputFile->Get("hNEventsProcessedPerRun");
2287 // round the number of events to avoid precision surprizes
2288 return TMath::Nint( hNEventsProcessedPerRun->GetBinContent(hNEventsProcessedPerRun->FindBin(run)) );
2291 //_________________________________________________________________________
2292 Long64_t GetTotalNumberOfEvents(Int_t nruns, Int_t runNumbers[])
2294 // Return total number of events in all runs
2296 Long64_t ntotal = 0;
2298 for (Int_t i = 0; i < nruns; i++)
2299 ntotal += GetNumberOfEvents(runNumbers[i]);
2304 //_________________________________________________________________________
2305 void DrawRunsDistribution(Int_t nruns, Int_t runNumbers[], Int_t dnbins = 100)
2307 // Draw events and events distribution;
2308 // dnbins -- number of bins in distribution.
2310 // initialize run index histogram ...
2311 TH1* hNEventsPerRunIndex = new TH1F("hNEventsPerRunIndex", "Number of processed events per run index", nruns,0,nruns);
2312 SetRunLabel(hNEventsPerRunIndex,nruns,runNumbers,1);
2313 // hNEventsPerRunIndex->SetXTitle("Run index");
2314 hNEventsPerRunIndex->SetYTitle("Number of events");
2317 for (Int_t i = 0; i < nruns; i++)
2318 hNEventsPerRunIndex->SetBinContent(i+1, GetNumberOfEvents(runNumbers[i]));
2320 // initialize distribution; 1.01 - to see the last bin
2321 TH1* distrib = new TH1F("hNEventsPerRunIndexDistr", "", dnbins, 0, hNEventsPerRunIndex->GetMaximum()*1.01);
2322 distrib->SetXTitle("Number of events");
2323 distrib->SetYTitle("Number of runs");
2325 // fill distribution
2326 for (Int_t i = 1; i <= nruns; i++)
2327 distrib->Fill(hNEventsPerRunIndex->GetBinContent(i));
2329 // draw histogram + distribution
2330 TCanvas *c1 = new TCanvas("hNEventsPerRunIndex","hNEventsPerRunIndex", 800,600);
2334 gPad->SetLeftMargin(0.06);
2335 gPad->SetRightMargin(0.04);
2336 gPad->SetTopMargin(0.10);
2337 gPad->SetBottomMargin(0.14);
2339 hNEventsPerRunIndex->SetTitleOffset(0.6,"Y");
2340 hNEventsPerRunIndex->SetTickLength(0.01,"Y");
2341 hNEventsPerRunIndex->SetLabelSize(0.06,"X");
2342 hNEventsPerRunIndex->Draw();
2345 gPad->SetLeftMargin(0.06);
2346 gPad->SetRightMargin(0.04);
2347 gPad->SetTopMargin(0.10);
2348 gPad->SetBottomMargin(0.13);
2349 distrib->SetTitleOffset(0.6,"Y");
2350 distrib->SetTickLength(0.01,"Y");
2356 //_________________________________________________________________________
2357 void ExcludeSmallRuns(Int_t &nruns, Int_t runNumbers[], Int_t nEventsMin = 100e+3)
2359 // Exclude runs with number of events < nEventsMin
2361 printf("Excluding runs (< %.0fk events):", nEventsMin/1000.);
2363 for (Int_t i = 0; i < nruns; i++)
2364 if (GetNumberOfEvents(runNumbers[i]) < nEventsMin) {
2365 printf(" %i", runNumbers[i]);
2367 runNumbers[i] = runNumbers[nruns];
2373 SortArray(nruns, runNumbers);
2376 //_________________________________________________________________________
2377 void ExcludeRunNumbers(Int_t &nruns, Int_t runNumbers[], Int_t nexcl, Int_t runs2Exclude[])
2379 // Exclude particular runs.
2381 // nexcl -- number of runs in runs2Exclude[];
2382 // runs2Exclude -- array with run numbers to exclude.
2384 for (Int_t i = 0; i < nruns; i++)
2385 for (Int_t e = 0; e < nexcl; e++)
2386 if (runNumbers[i] == runs2Exclude[e]) {
2387 Printf("Excluding run: %i", runs2Exclude[e]);
2389 runNumbers[i] = runNumbers[nruns];
2394 SortArray(nruns, runNumbers);
2397 //_________________________________________________________________________
2398 void SortArray(const Int_t nruns, Int_t runNumbers[])
2401 // used for runNumbers array sorting.
2403 Int_t indices[nruns];
2404 Int_t runNumbers_unsort[nruns];
2406 for (Int_t i = 0; i < nruns; i++)
2407 runNumbers_unsort[i] = runNumbers[i];
2409 TMath::Sort(nruns, runNumbers_unsort, indices, kFALSE);
2411 for (Int_t i = 0; i < nruns; i++)
2412 runNumbers[i] = runNumbers_unsort[indices[i]];
2415 //_________________________________________________________________________
2416 void PrintRunNumbers(Int_t nruns, Int_t runNumbers[])
2418 // Print a table run index / run number / number of events.
2421 Printf("| index | run number | nevents |");
2422 Printf("|-------|------------|-----------|");
2424 for (Int_t i = 0; i < nruns; i++)
2425 Printf("| %-4i | %-9i | %-8i |", i, runNumbers[i], GetNumberOfEvents(runNumbers[i]));
2427 Printf("| Events in total: %-10lli |\n", GetTotalNumberOfEvents(nruns, runNumbers));
2430 //------------------------------------------------------------------------
2431 void SetRunLabel(TH1 *histo, Int_t nruns, Int_t *runNumbers, Int_t axis)
2434 for (Int_t i=0; i<nruns; i++) {
2435 runText = Form("%d",runNumbers[i]);
2437 histo->GetXaxis()->SetBinLabel(i+1,runText.Data());
2439 histo->GetYaxis()->SetBinLabel(i+1,runText.Data());
2441 histo->LabelsOption("v");