Belle II Software development
DQMHistAnalysisECLSummary.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9//THIS MODULE
10#include <dqm/analysis/modules/DQMHistAnalysisECLSummary.h>
11
12//ECL
13#include <ecl/geometry/ECLGeometryPar.h>
14
15//ROOT
16#include <TStyle.h>
17#include <TLine.h>
18
19//boost
20#include "boost/format.hpp"
21
22using namespace Belle2;
23
24REG_MODULE(DQMHistAnalysisECLSummary);
25
28 m_neighbours_obj("F", 0.1)
29{
30
31 B2DEBUG(20, "DQMHistAnalysisECLSummary: Constructor done.");
32 addParam("pvPrefix", m_pvPrefix, "Prefix to use for PVs registered by this module",
33 std::string("ECL:channels_info:"));
34 addParam("useChannelMask", m_useChannelMask,
35 "Mask Cell IDs based on information from ECL PVs",
36 true);
37 addParam("maxDeviationForOccupancy", m_maxDeviationForOccupancy,
38 "The higher this parameter, the larger differences in occupancy are allowed for adjacent channels", 0.28);
39 addParam("maxDeviationForChi2", m_maxDeviationForChi2,
40 "The higher this parameter, the larger differences in the number of hits with bad chi2 are allowed for adjacent channels",
41 2.5);
42 addParam("onlyIfUpdated", m_onlyIfUpdated, "If true (default), update EPICS PVs only if histograms were updated.",
43 true);
44}
45
46
50
52{
53 m_mapper.initFromFile();
54
55 //=== Set up ECL alarms and corresponding PVs
56
57 m_ecl_alarms = {
58 {"dead", "#splitline{dead}{channels}", 1, 1, 1e5},
59 {"cold", "#splitline{cold}{channels}", 1, 2, 1e5},
60 {"hot", "#splitline{hot}{channels}", 25, 50, 1e5},
61 {"bad_chi2", "#splitline{bad #chi^{2}}{channels}", 5, 10, 1e6},
62 {"bad_fit", "#splitline{fit incon-}{sistencies}", 5, 10, 0 },
63 };
64
65 // Prepare EPICS PVs
66 for (auto& alarm : m_ecl_alarms) {
67 // By crate
68 for (int crate_id = 1; crate_id <= ECL::ECL_CRATES; crate_id++) {
69 std::string pv_name = (boost::format("crate%02d:%s") % crate_id % alarm.name).str();
70 registerEpicsPV(m_pvPrefix + pv_name, pv_name);
71 }
72 // Totals
73 for (auto& ecl_part : {"All", "FWDEndcap", "Barrel", "BWDEndcap"}) {
74 std::string pv_name = (boost::format("%s:%s") % ecl_part % alarm.name).str();
75 registerEpicsPV(m_pvPrefix + pv_name, pv_name);
76 }
77 // Masked Cell IDs
78 std::string mask_pv_name = (boost::format("mask:%s") % alarm.name).str();
79 registerEpicsPV(m_pvPrefix + mask_pv_name, mask_pv_name);
80 }
81
83
84 //=== Set up the histogram to indicate alarm status
85
86 TString title = "#splitline{ECL errors monitoring}";
87 title += "{E - Error, W - Warning, L - Low statistics}";
88 title += ";ECLCollector ID (same as Crate ID)";
89 h_channels_summary = new TH2F("channels_summary", title,
90 ECL::ECL_CRATES, 1, ECL::ECL_CRATES + 1,
91 m_ecl_alarms.size(), 0, m_ecl_alarms.size());
92
93 // Do not show statistics box.
94 h_channels_summary->SetStats(0);
95 h_channels_summary->SetMinimum(0);
96 h_channels_summary->SetMaximum(1);
97
98 //=== Set X axis labels
99
100 for (int i = 1; i <= ECL::ECL_CRATES; i++) {
101 h_channels_summary->GetXaxis()->SetBinLabel(i, std::to_string(i).c_str());
102 }
103 h_channels_summary->LabelsOption("v", "X"); // Rotate X axis labels 90 degrees
104 h_channels_summary->SetTickLength(0, "XY");
105
106 //=== Customize offsets and margins
107
108 h_channels_summary->GetXaxis()->SetTitleOffset(0.95);
109 h_channels_summary->GetXaxis()->SetTitleSize(0.05);
110 h_channels_summary->GetXaxis()->SetLabelSize(0.04);
111 h_channels_summary->GetYaxis()->SetLabelSize(0.06);
112
113 c_channels_summary = new TCanvas("ECL/c_channels_summary_analysis");
114 c_channels_summary->SetTopMargin(0.10);
115 c_channels_summary->SetLeftMargin(0.20);
116 c_channels_summary->SetRightMargin(0.005);
117 c_channels_summary->SetBottomMargin(0.10);
118
119 //=== Additional canvases/histograms to display which channels have problems
120 c_occupancy = new TCanvas("ECL/c_cid_Thr5MeV_overlaid_analysis");
121 c_bad_chi2 = new TCanvas("ECL/c_bad_quality_overlaid_analysis");
122
123 h_bad_occ_overlay = new TH1F("bad_occ_overlay", "", ECL::ECL_TOTAL_CHANNELS,
124 1, ECL::ECL_TOTAL_CHANNELS + 1);
125 h_bad_occ_overlay->SetLineColor(kRed);
126 h_bad_occ_overlay->SetLineStyle(kDashed);
127 h_bad_occ_overlay->SetFillColor(kRed);
128 h_bad_occ_overlay->SetFillStyle(3007);
129
130 h_bad_chi2_overlay = new TH1F("bad_chi2_overlay", "", ECL::ECL_TOTAL_CHANNELS,
131 1, ECL::ECL_TOTAL_CHANNELS + 1);
132 h_bad_chi2_overlay->SetLineColor(kRed);
133 h_bad_chi2_overlay->SetLineStyle(kDashed);
134 h_bad_chi2_overlay->SetFillColor(kRed);
135 h_bad_chi2_overlay->SetFillStyle(3007);
136
137 h_bad_occ_overlay_green = new TH1F("bad_occ_overlay_green", "", ECL::ECL_TOTAL_CHANNELS,
138 1, ECL::ECL_TOTAL_CHANNELS + 1);
139 h_bad_occ_overlay_green->SetLineColor(kGreen);
140 h_bad_occ_overlay_green->SetLineStyle(kDotted);
141 h_bad_occ_overlay_green->SetFillColor(kGreen);
142 h_bad_occ_overlay_green->SetFillStyle(3013);
143
144 h_bad_chi2_overlay_green = new TH1F("bad_chi2_overlay_green", "", ECL::ECL_TOTAL_CHANNELS,
145 1, ECL::ECL_TOTAL_CHANNELS + 1);
146 h_bad_chi2_overlay_green->SetLineColor(kGreen);
147 h_bad_chi2_overlay_green->SetLineStyle(kDotted);
148 h_bad_chi2_overlay_green->SetFillColor(kGreen);
149 h_bad_chi2_overlay_green->SetFillStyle(3013);
150
151 B2DEBUG(20, "DQMHistAnalysisECLSummary: initialized.");
152}
153
155{
156 B2DEBUG(20, "DQMHistAnalysisECLSummary: beginRun called.");
157
158 //=== Update m_ecl_alarms based on PV limits
159
161
162 //=== Set Y axis labels from m_ecl_alarms
163
164 for (size_t i = 0; i < m_ecl_alarms.size(); i++) {
165 TString label = m_ecl_alarms[i].title;
166 label += " #geq ";
167 label += m_ecl_alarms[i].alarm_limit;
168 h_channels_summary->GetYaxis()->SetBinLabel(i + 1, label);
169 }
170}
171
173{
174 TH1* h_total_events = findHist("ECL/event", m_onlyIfUpdated);
175 if (!h_total_events) return;
176 m_total_events = h_total_events->GetEntries();
177
178 c_channels_summary->Clear();
179 h_channels_summary->Reset();
180
181 // [alarm_type][crate_id - 1] -> number of problematic channels in that crate
182 std::vector< std::vector<int> > alarm_counts = updateAlarmCounts();
183
184 //=== Set warning and error indicators on the histogram
185
186 const double HISTCOLOR_RED = 0.9;
187 const double HISTCOLOR_GREEN = 0.45;
188 const double HISTCOLOR_ORANGE = 0.65;
189 const double HISTCOLOR_GRAY = 0.00;
190
191 // TODO: Currently, I am explicitly removing
192 // TText labels from the histogram. However,
193 // normally I would assume that canvas->Clear()
194 // should have already removed them. So
195 // it might indicate a memory leak or some
196 // problem with ROOT->JSON->ROOT conversion
197 for (auto& text : m_labels) {
198 delete text;
199 }
200 m_labels.clear();
201
202 //== Set correct warning/error colors for each bin based on alarm thresholds
203
204 bool enough = false;
205 bool warn_flag = false;
206 bool error_flag = false;
207
208 for (size_t alarm_idx = 0; alarm_idx < alarm_counts.size(); alarm_idx++) {
209 for (size_t crate = 0; crate < alarm_counts[alarm_idx].size(); crate++) {
210 double color;
211 char label_text[2] = {0};
212 int alarm_limit = m_ecl_alarms[alarm_idx].alarm_limit;
213 int warning_limit = m_ecl_alarms[alarm_idx].warning_limit;
214
215 int alarms = alarm_counts[alarm_idx][crate];
216
217 if (m_total_events < m_ecl_alarms[alarm_idx].required_statistics || alarms < 0) {
218 // `alarms < 0` means that we could not get the relevant histogram.
219 // We assume that the histogram is not yet available due to low
220 // statistics but this is not a guarantee.
221 color = HISTCOLOR_GRAY;
222 label_text[0] = 'L';
223 } else if (alarms >= alarm_limit) {
224 color = HISTCOLOR_RED;
225 label_text[0] = 'E';
226 error_flag = true;
227 enough = true;
228 } else if (alarms >= warning_limit) {
229 color = HISTCOLOR_ORANGE;
230 label_text[0] = 'W';
231 warn_flag = true;
232 enough = true;
233 } else {
234 color = HISTCOLOR_GREEN;
235 enough = true;
236 }
237 if (label_text[0] == 'E' || label_text[0] == 'W') {
238 B2DEBUG(100, "Non-zero (" << alarm_counts[alarm_idx][crate]
239 << ") for alarm_idx, crate = " << alarm_idx << ", " << crate);
240 }
241 h_channels_summary->SetBinContent(crate + 1, alarm_idx + 1, color);
242 if (label_text[0] != 0) {
243 auto text = new TText((crate + 1.5), (alarm_idx + 0.5), label_text);
244 if (label_text[0] != 'L') text->SetTextColor(kWhite);
245 text->SetTextSize(0.03);
246 text->SetTextAlign(22); // centered
247 m_labels.push_back(text);
248 }
249 }
250 }
251
252 //=== Draw histogram, labels and grid
253
254 // Customize title
255 auto gstyle_title_h = gStyle->GetTitleH();
256 auto gstyle_title_x = gStyle->GetTitleX();
257 auto gstyle_title_y = gStyle->GetTitleY();
258 gStyle->SetTitleH(0.04);
259 gStyle->SetTitleX(0.60);
260 gStyle->SetTitleY(1.00);
261
262 //=== Set background color based on combined status.
263 c_channels_summary->cd();
264 colorizeCanvas(c_channels_summary, makeStatus(enough, warn_flag, error_flag));
266
267 //=== Prepare special style objects to use correct color palette
268 // and use it only for this histogram
269 // Based on https://root-forum.cern.ch/t/different-color-palettes-for-different-plots-with-texec/5250/3
270 // and https://root.cern/doc/master/multipalette_8C.html
271
272 //=== Draw with special style
273 // https://root.cern.ch/js/latest/examples.htm#th2_colpal77
274 h_channels_summary->Draw("");
275 h_channels_summary->Draw("col;pal55;same");
276 for (auto& text : m_labels) {
277 text->Draw();
278 }
280
281 //
282 c_channels_summary->Modified();
283 c_channels_summary->Update();
284 c_channels_summary->Draw();
285
286 m_ecl_style = new TExec("ecl_style",
287 "gStyle->SetPalette(kRainBow);");
288 c_channels_summary->GetListOfPrimitives()->AddFirst(m_ecl_style);
289 m_default_style = new TExec("default_style",
290 "gStyle->SetPalette(kBird);"); // " Changing back to the default color palette
291 c_channels_summary->GetListOfPrimitives()->AddLast(m_default_style);
292
293 gStyle->SetTitleH(gstyle_title_h);
294 gStyle->SetTitleX(gstyle_title_x);
295 gStyle->SetTitleY(gstyle_title_y);
296}
297
299{
300 B2DEBUG(20, "DQMHistAnalysisECLSummary: endRun called");
301
302 updateAlarmCounts(true);
303}
304
305
307{
308 B2DEBUG(20, "terminate called");
309 delete c_channels_summary;
310 delete c_occupancy;
311 delete c_bad_chi2;
312 delete h_bad_occ_overlay;
313 delete h_bad_chi2_overlay;
314}
315
316std::pair<int, DQMHistAnalysisECLSummaryModule::ECLAlarmType> DQMHistAnalysisECLSummaryModule::getAlarmByName(std::string name)
317{
318 int index = 0;
319 for (auto& alarm_info : m_ecl_alarms) {
320 if (alarm_info.name == name) return {index, alarm_info};
321 index++;
322 }
323 B2FATAL("Could not get ECL alarm " + name);
324 return {-1, m_ecl_alarms[0]};
325}
326
328{
329 //===== Get alarm limits
330
331 for (auto& alarm : m_ecl_alarms) {
332 // In the current version, use only first crate PV to get alarm limits
333 int crate_id = 1;
334
335 std::string pv_name = (boost::format("crate%02d:%s") % crate_id % alarm.name).str();
336 double unused = 0;
337 double upper_warning_limit = 0;
338 double upper_alarm_limit = 0;
339 bool accessed = requestLimitsFromEpicsPVs(pv_name, unused, unused, upper_warning_limit, upper_alarm_limit);
340
341 if (!accessed || upper_alarm_limit <= 0 || upper_warning_limit <= 0) {
342 B2WARNING("Failed to get alarm limits");
343 return;
344 }
345
346 alarm.alarm_limit = upper_alarm_limit;
347 alarm.warning_limit = upper_warning_limit;
348 }
349
350 //===== Get masked channels
351
352 // alarm_type -> array_of_masked_channels
353 static std::map<std::string, dbr_sts_long_array> mask_info;
354 if (!getMaskedChannels(mask_info)) {
355 B2WARNING("Failed to get arrays of masked channels");
356 return;
357 }
358
359 for (auto& alarm : m_ecl_alarms) {
360 // Update the list of masked channels
361 m_mask[alarm.name].clear();
362 for (int i = 0; i < c_max_masked_channels; i++) {
363 int cell_id = mask_info[alarm.name].value[i];
364 // The array of masked channels goes like this: [15, 22, 716, 0, 0, 0, 0, 0]
365 // So if 0 (invalid channel ID) is encountered, stop reading further.
366 if (cell_id == 0) break;
367 m_mask[alarm.name].insert(cell_id);
368 }
369 }
370}
371
372bool DQMHistAnalysisECLSummaryModule::getMaskedChannels(std::map<std::string, dbr_sts_long_array>& mask_info)
373{
374 for (auto& alarm : m_ecl_alarms) {
375 std::string mask_pv_name = (boost::format("mask:%s") % alarm.name).str();
376 chid mask_chid = getEpicsPVChID(mask_pv_name);
377
378 if (mask_chid == nullptr) return false;
379 // ca_array_get data is NOT VALID here.
380 // It will only be valid after ca_pend_io returns ECA_NORMAL
381 auto r = ca_array_get(DBR_STS_LONG, c_max_masked_channels, mask_chid, &mask_info[alarm.name]);
382 if (r != ECA_NORMAL) return false;
383 }
384
385 // ca_pend_io has to be called here,
386 // see https://epics.anl.gov/base/R3-14/12-docs/CAref.html#ca_get
387
388 auto r = ca_pend_io(5.0);
389 if (r != ECA_NORMAL) return false;
390
391 return true;
392}
393
394std::vector< std::vector<int> > DQMHistAnalysisECLSummaryModule::updateAlarmCounts(bool update_mirabelle)
395{
396 std::vector< std::vector<int> > alarm_counts(
397 m_ecl_alarms.size(), std::vector<int>(ECL::ECL_CRATES));
398
399 //=== Get number of fit inconsistencies
400
401 TH1* h_fail_crateid = findHist("ECL/fail_crateid", false);
402
403 const int fit_alarm_index = getAlarmByName("bad_fit").first;
404 for (int crate_id = 1; crate_id <= ECL::ECL_CRATES; crate_id++) {
405 int errors_count = 0;
406 if (h_fail_crateid) {
407 errors_count = h_fail_crateid->GetBinContent(crate_id);
408 } else {
409 errors_count = -1;
410 }
411
412 alarm_counts[fit_alarm_index][crate_id - 1] = errors_count;
413 }
414
415 // [Cell ID] -> error_bitmask
416 std::map<int, int> error_bitmasks;
417
418 //=== Get number of dead/cold/hot channels
419 // cppcheck-suppress unassignedVariable
420 for (auto& [cell_id, error_bitmask] : getChannelsWithOccupancyProblems()) {
421 error_bitmasks[cell_id] |= error_bitmask;
422 }
423 //=== Get number of channels with bad_chi2
424 // cppcheck-suppress unassignedVariable
425 for (auto& [cell_id, error_bitmask] : getChannelsWithChi2Problems()) {
426 error_bitmasks[cell_id] |= error_bitmask;
427 }
428
429 //=== Combine the information
430
431 if (!update_mirabelle) {
432 h_bad_chi2_overlay->Reset();
433 h_bad_occ_overlay->Reset();
436 }
437
438 static std::vector<std::string> indices = {"dead", "cold", "hot", "bad_chi2"};
439 for (auto& index_name : indices) {
440 int alarm_index = getAlarmByName(index_name).first;
441 int alarm_bit = 1 << alarm_index;
442
443 for (auto& [cid, error_bitmask] : error_bitmasks) {
444 if ((error_bitmask & alarm_bit) == 0) continue;
445
446 bool masked = (m_mask[index_name].find(cid) != m_mask[index_name].end());
447
448 if (!masked) {
449 int crate_id = m_mapper.getCrateID(cid);
450 // If channel is not masked, increase the alarm counter.
451 alarm_counts[alarm_index][crate_id - 1] += 1;
452 }
453
454 if (update_mirabelle) continue;
455
456 // Unless it is the end of run, add problematic channels
457 // to the detailed 1D histogram.
458
459 if (index_name == "bad_chi2") {
460 if (!masked) h_bad_chi2_overlay->SetBinContent(cid, 1);
461 else h_bad_chi2_overlay_green->SetBinContent(cid, 1);
462 } else if (index_name == "dead" || index_name == "cold" || index_name == "hot") {
463 if (!masked) h_bad_occ_overlay->SetBinContent(cid, 1);
464 else h_bad_occ_overlay_green->SetBinContent(cid, 1);
465 }
466 }
467
468 if (update_mirabelle) continue;
469
470 //=== Prepare the histogram to display the list of all bad channels
471 if (index_name == "hot" || index_name == "bad_chi2") {
472 TCanvas* current_canvas;
473 TH1* main_hist;
474 TH1F* overlay_hist;
475 TH1F* overlay_hist_green;
476 if (index_name == "hot") {
477 main_hist = findHist("ECL/cid_Thr5MeV", m_onlyIfUpdated);
478 overlay_hist = h_bad_occ_overlay;
479 overlay_hist_green = h_bad_occ_overlay_green;
480 current_canvas = c_occupancy;
481 } else {
482 main_hist = findHist("ECL/bad_quality", m_onlyIfUpdated);
483 overlay_hist = h_bad_chi2_overlay;
484 overlay_hist_green = h_bad_chi2_overlay_green;
485 current_canvas = c_bad_chi2;
486 }
487
488 if (main_hist && main_hist->GetEntries() > 0) {
489 for (auto& overlay : {overlay_hist, overlay_hist_green}) {
490 for (int bin_id = 1; bin_id <= ECL::ECL_TOTAL_CHANNELS; bin_id++) {
491 if (overlay->GetBinContent(bin_id) == 0) continue;
492 // Do not adjust bin height for dead channels
493 if (main_hist->GetBinContent(bin_id) < 1e-6) continue;
494 overlay->SetBinContent(bin_id, main_hist->GetBinContent(bin_id));
495 }
496 }
497
498 current_canvas->Clear();
499 current_canvas->cd();
500 main_hist->Draw("hist");
501 overlay_hist->Draw("hist;same");
502 overlay_hist_green->Draw("hist;same");
503 current_canvas->Modified();
504 current_canvas->Update();
505 current_canvas->Draw();
506 }
507 }
508 }
509
510 //== Update EPICS PVs or MiraBelle monObjs
511
512 for (size_t alarm_idx = 0; alarm_idx < alarm_counts.size(); alarm_idx++) {
513 auto& alarm = m_ecl_alarms[alarm_idx];
514 std::map<std::string, int> total;
515 // Convert values per crate to totals
516 for (size_t crate = 0; crate < alarm_counts[alarm_idx].size(); crate++) {
517 int crate_id = crate + 1;
518 int value = alarm_counts[alarm_idx][crate];
519
520 if (!update_mirabelle) {
521 std::string pv_name = (boost::format("crate%02d:%s") % crate_id % alarm.name).str();
522 setEpicsPV(pv_name, value);
523 }
524
525 total["All"] += value;
526 if (crate_id <= ECL::ECL_BARREL_CRATES) {
527 total["Barrel"] += value;
528 } else if (crate_id <= ECL::ECL_BARREL_CRATES + ECL::ECL_FWD_CRATES) {
529 total["FWDEndcap"] += value;
530 } else {
531 total["BWDEndcap"] += value;
532 }
533 }
534 // Export totals
535 for (auto& ecl_part : {"All", "FWDEndcap", "Barrel", "BWDEndcap"}) {
536 std::string pv_name = (boost::format("%s:%s") % ecl_part % alarm.name).str();
537 if (update_mirabelle) {
538 std::string var_name = pv_name;
539 std::replace(var_name.begin(), var_name.end(), ':', '_');
540 B2DEBUG(100, var_name << " = " << total[ecl_part]);
541 m_monObj->setVariable(var_name, total[ecl_part]);
542 } else {
543 setEpicsPV(pv_name, total[ecl_part]);
544 }
545 }
546 }
547
548 return alarm_counts;
549}
550
552{
553 static std::vector< std::vector<short> > neighbours(ECL::ECL_TOTAL_CHANNELS);
554 if (neighbours[0].size() == 0) {
555 for (int cid0 = 0; cid0 < ECL::ECL_TOTAL_CHANNELS; cid0++) {
556 // [0]First is the crystal itself.
557 // [1]Second is PHI neighbour, phi+1.
558 // [2]Third is PHI neighbour, phi-1.
559 // [3]Next one (sometimes two) are THETA neighbours in theta-1.
560 // [4]Next one (sometimes two) are THETA neighbours in theta+1.
561 neighbours[cid0] = m_neighbours_obj.getNeighbours(cid0 + 1);
562 // Remove first element (the crystal itself)
563 neighbours[cid0].erase(neighbours[cid0].begin());
564 }
565 }
566
567 TH1* h_occupancy = findHist("ECL/cid_Thr5MeV", m_onlyIfUpdated);
568 const double max_deviation = m_maxDeviationForOccupancy;
569 return getSuspiciousChannels(h_occupancy, m_total_events, neighbours,
570 max_deviation, true);
571}
572
574{
576
577 static std::vector< std::vector<short> > neighbours(ECL::ECL_TOTAL_CHANNELS);
578 if (neighbours[0].size() == 0) {
579 for (int cid_center = 1; cid_center <= ECL::ECL_TOTAL_CHANNELS; cid_center++) {
580 geom->Mapping(cid_center - 1);
581 const int theta_id_center = geom->GetThetaID();
582 int phi_id_center = geom->GetPhiID();
583 const int crystals_in_ring = m_neighbours_obj.getCrystalsPerRing(theta_id_center);
584 phi_id_center = phi_id_center * 144 / crystals_in_ring;
585 for (int cid0 = 0; cid0 < 8736; cid0++) {
586 if (cid0 == cid_center - 1) continue;
587 geom->Mapping(cid0);
588 int theta_id = geom->GetThetaID();
589 int phi_id = geom->GetPhiID();
590 phi_id = phi_id * 144 / m_neighbours_obj.getCrystalsPerRing(theta_id);
591 if (std::abs(theta_id - theta_id_center) <= 2 &&
592 std::abs(phi_id - phi_id_center) <= 2) {
593 neighbours[cid_center - 1].push_back(cid0);
594 }
595 }
596 }
597 }
598
599 TH1* h_bad_chi2 = findHist("ECL/bad_quality", m_onlyIfUpdated);
600 const double max_deviation = m_maxDeviationForChi2;
601 return getSuspiciousChannels(h_bad_chi2, m_total_events, neighbours,
602 max_deviation, false);
603}
604
606 TH1* hist, double total_events,
607 const std::vector< std::vector<short> >& neighbours,
608 double max_deviation, bool occupancy_histogram)
609{
610 std::map<int, int> retval;
611
612 if (!hist) return retval;
613 //== A bit hacky solution to skip incorrectly
614 // filled histograms
615 if (hist->Integral() <= 0) return retval;
616
617 //=== Extract alarm details
618 // cppcheck-suppress unassignedVariable
619 const auto& [dead_index, dead_alarm] = getAlarmByName("dead");
620 // cppcheck-suppress unassignedVariable
621 const auto& [cold_index, cold_alarm] = getAlarmByName("cold");
622 // cppcheck-suppress unassignedVariable
623 const auto& [hot_index, hot_alarm ] = getAlarmByName("hot");
624 // cppcheck-suppress unassignedVariable
625 const auto& [chi2_index, chi2_alarm] = getAlarmByName("bad_chi2");
626
627 double min_required_events;
628
629 if (occupancy_histogram) {
630 min_required_events = std::min({
631 dead_alarm.required_statistics,
632 cold_alarm.required_statistics,
633 hot_alarm.required_statistics,
634 });
635 } else {
636 min_required_events = chi2_alarm.required_statistics;
637 }
638
639 if (total_events < min_required_events) return retval;
640
641 int dead_bit = 1 << dead_index;
642 int cold_bit = 1 << cold_index;
643 int hot_bit = 1 << hot_index;
644 int chi2_bit = 1 << chi2_index;
645
646 // == Search for dead channels
647
648 if (occupancy_histogram) {
649 // This indicates that DQMHistAnalysisECL module is not included in the path
650 bool not_normalized = (findCanvas("ECL/c_cid_Thr5MeV_analysis") == nullptr);
651 if (total_events >= dead_alarm.required_statistics) {
652 double min_occupancy;
653 const std::string run_type = getRunType();
654 if (run_type == "null" || run_type == "debug" || run_type == "cosmic") {
655 // For null runs, occupancy should be higher than 0.0001%
656 min_occupancy = 1e-6;
657 } else if (run_type == "physics") {
658 // For physics runs, occupancy should be higher than 0.01%
659 min_occupancy = 1e-4;
660 }
661 if (not_normalized) {
662 // The histogram is not normalized, multiply the threshold by evt count
663 min_occupancy *= total_events;
664 }
665 for (int cid = 1; cid <= ECL::ECL_TOTAL_CHANNELS; cid++) {
666 if (hist->GetBinContent(cid) > min_occupancy) continue;
667 retval[cid] |= dead_bit;
668 }
669 }
670 if (total_events >= hot_alarm.required_statistics) {
671 // Number of hits with E > threshold should be less than 70%
672 double max_occupancy = 0.7;
673 if (not_normalized) {
674 // The histogram is not normalized, multiply the threshold by evt count
675 max_occupancy *= total_events;
676 }
677 for (int cid = 1; cid <= ECL::ECL_TOTAL_CHANNELS; cid++) {
678 if (hist->GetBinContent(cid) < max_occupancy) continue;
679 retval[cid] |= hot_bit;
680 }
681 }
682 }
683
684 // == Search for cold and hot channels (or channels with bad chi2)
685
686 for (int cid = 1; cid <= ECL::ECL_TOTAL_CHANNELS; cid++) {
687 double actual_value = hist->GetBinContent(cid);
688
689 std::vector<short> neighb = neighbours[cid - 1];
690 std::multiset<double> values_sorted;
691 for (auto& neighbour_cid : neighb) {
692 values_sorted.insert(hist->GetBinContent(neighbour_cid));
693 }
694
695 double expected_value;
696 if (occupancy_histogram) {
697 // Use "lower 25%" value:
698 expected_value = *std::next(values_sorted.begin(), 1 * neighb.size() / 4);
699 } else {
700 // Use "upper ~75%" value (48 is the expected number of neighbours)
701 expected_value = *std::next(values_sorted.begin(), 35 * neighb.size() / 48);
702 }
703 double deviation = std::abs((actual_value - expected_value) /
704 expected_value);
705
706 if (!occupancy_histogram) {
707 // Min occupancy for high-energy (> 1GeV) hits
708 double min_occupancy = 1.41e-5;
709 if (findCanvas("ECL/c_bad_quality_analysis") == nullptr) {
710 // The histogram is not normalized, multiply the threshold by evt count
711 min_occupancy *= total_events;
712 }
713
714 if (actual_value < min_occupancy) continue;
715 }
716
717 if (deviation < max_deviation) continue;
718
719 if (occupancy_histogram) {
720 if (actual_value < expected_value) retval[cid] |= cold_bit;
721 if (actual_value > 2 * expected_value) retval[cid] |= hot_bit;
722 } else {
723 if (actual_value > expected_value) retval[cid] |= chi2_bit;
724 }
725 }
726
727 return retval;
728}
729
731{
732 static std::map<std::string, std::vector<TLine*> > lines;
733 std::string name = hist->GetName();
734 if (lines[name].empty()) {
735 int x_min = hist->GetXaxis()->GetXmin();
736 int x_max = hist->GetXaxis()->GetXmax();
737 int y_min = hist->GetYaxis()->GetXmin();
738 int y_max = hist->GetYaxis()->GetXmax();
739 for (int x = x_min + 1; x < x_max; x++) {
740 auto l = new TLine(x, 0, x, 5);
741 l->SetLineStyle(kDashed);
742 lines[name].push_back(l);
743 }
744 for (int y = y_min + 1; y < y_max; y++) {
745 auto l = new TLine(1, y, ECL::ECL_CRATES + 1, y);
746 l->SetLineStyle(kDashed);
747 lines[name].push_back(l);
748 }
749 }
750 for (auto line : lines[name]) {
751 line->Draw();
752 }
753}
754
std::map< int, int > getSuspiciousChannels(TH1 *hist, double total_events, const std::vector< std::vector< short > > &neighbours, double max_deviation, bool occupancy_histogram)
Get outlier channels, ones with values (occupancy, errors, etc) that are much higher than the values ...
double m_total_events
Number of events with ECL data in the current run.
double m_maxDeviationForChi2
The higher this parameter, the larger differences in the number of hits with bad chi2 are allowed for...
TH1F * h_bad_occ_overlay
Overlay to indicate bad channels on occupancy histogram.
void initialize() override final
Initialize the module.
void updateAlarmConfig()
Set alarm limits in DQM based on EPICS PV limits.
std::vector< std::vector< int > > updateAlarmCounts(bool update_mirabelle=false)
Calculate and return number of problems per crate for all alarm categories.
std::vector< TText * > m_labels
Labels on the ECL alarms histogram.
ECL::ECLChannelMapper m_mapper
Object to map ECL crystal ID to ECL crate ID.
void drawGrid(TH2 *hist)
Draw grid with TLine primitives for the specified histogram.
TCanvas * c_channels_summary
TCanvas for ECL alarms regarding suspicious channels.
static const int c_max_masked_channels
Size of an array with masked channels.
std::pair< int, ECLAlarmType > getAlarmByName(std::string name)
Returns index of the specific alarm type and detailed information.
TH2F * h_channels_summary
Summarized information about ECL channels.
double m_maxDeviationForOccupancy
The higher this parameter, the larger differences in occupancy are allowed for adjacent channels.
ECL::ECLNeighbours m_neighbours_obj
Object to get ECL crystal neighbours.
TCanvas * c_occupancy
ECL occupancy histogram with highlighted suspicious channels.
TExec * m_ecl_style
Special object to specify style parameters.
std::string m_pvPrefix
Prefix to use for PVs registered by this module.
MonitoringObject * m_monObj
MiraBelle monitoring object.
void event() override final
Event processor.
TCanvas * c_bad_chi2
ECL bad chi2 histogram with highlighted suspicious channels.
TH1F * h_bad_chi2_overlay
Overlay to indicate bad channels on chi2 histogram.
void endRun() override final
Call when a run ends.
bool m_onlyIfUpdated
If true (default), update EPICS PVs only if there were changes in the histograms.
TH1F * h_bad_chi2_overlay_green
Overlay to indicate masked bad channels on chi2 histogram.
void beginRun() override final
Call when a run begins.
DQMHistAnalysisECLSummaryModule()
< derived from DQMHistAnalysisModule class.
TExec * m_default_style
Special object to revert changes done by ecl_style.
TH1F * h_bad_occ_overlay_green
Overlay to indicate masked bad channels on occupancy histogram.
std::map< std::string, std::set< int > > m_mask
List of masked Cell IDs for each alarm type.
bool getMaskedChannels(std::map< std::string, dbr_sts_long_array > &mask_info)
Get the array of masked channels for each type of alarm case monitored by this module.
bool m_useChannelMask
If true, mask Cell IDs based on DQM:ECL:channels_info:{alarm_type} PV info.
std::vector< ECLAlarmType > m_ecl_alarms
Different alarms monitored in h_channels_summary.
TCanvas * findCanvas(TString cname)
Find canvas by name.
static MonitoringObject * getMonitoringObject(const std::string &name)
Get MonitoringObject with given name (new object is created if non-existing)
void colorizeCanvas(TCanvas *canvas, EStatus status)
Helper function for Canvas colorization.
static TH1 * findHist(const std::string &histname, bool onlyIfUpdated=false)
Get histogram from list (no other search).
@ c_ColorTooFew
Not enough entries/event to judge.
static const std::string & getRunType(void)
Get the list of the reference histograms.
void setEpicsPV(std::string keyname, double value)
Write value to a EPICS PV.
DQMHistAnalysisModule()
Constructor / Destructor.
EStatus makeStatus(bool enough, bool warn_flag, bool error_flag)
Helper function to judge the status for coloring and EPICS.
int registerEpicsPV(std::string pvname, std::string keyname="")
EPICS related Functions.
chid getEpicsPVChID(std::string keyname)
Get EPICS PV Channel Id.
bool requestLimitsFromEpicsPVs(chid id, double &lowerAlarm, double &lowerWarn, double &upperWarn, double &upperAlarm)
Get Alarm Limits from EPICS PV.
The Class for ECL Geometry Parameters.
static ECLGeometryPar * Instance()
Static method to get a reference to the ECLGeometryPar instance.
void addParam(const std::string &name, T &paramVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
Definition Module.h:559
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition Module.h:649
Abstract base class for different kinds of events.