16 #include <ecl/modules/eclDQM/eclDQM.h>
19 #include <boost/format.hpp>
20 #include <boost/range/combine.hpp>
23 #include <framework/core/HistoModule.h>
24 #include <framework/dataobjects/EventMetaData.h>
25 #include <framework/gearbox/Unit.h>
26 #include <framework/logging/Logger.h>
29 #include <ecl/dataobjects/ECLDigit.h>
30 #include <ecl/dataobjects/ECLCalDigit.h>
31 #include <ecl/dataobjects/ECLTrig.h>
32 #include <ecl/dataobjects/ECLDsp.h>
33 #include <ecl/utility/ECLChannelMapper.h>
34 #include <ecl/geometry/ECLGeometryPar.h>
37 #include <mdst/dataobjects/TRGSummary.h>
43 #include <TDirectory.h>
57 m_calibrationThrApsd("ECL_FPGA_StoreWaveform")
60 setDescription(
"ECL Data Quality Monitor");
61 setPropertyFlags(c_ParallelProcessingCertified);
63 m_WaveformOption = {
"psd",
"logic",
"rand",
"dphy",
"other"};
65 addParam(
"histogramDirectoryName", m_histogramDirectoryName,
66 "histogram directory in ROOT file", std::string(
"ECL"));
67 addParam(
"EnergyUpperThr", m_EnergyUpperThr,
"Upper threshold of energy deposition in event, [GeV]", 20.0 *
Belle2::Unit::GeV);
68 addParam(
"HitThresholds", m_HitThresholds,
"Thresholds to display hit occupancy, MeV", std::vector<double> {0, 5, 10, 50});
69 addParam(
"TotalEnergyThresholds", m_TotalEnergyThresholds,
"Thresholds to display total energy, MeV", std::vector<double> {0, 5, 7});
70 addParam(
"TimingThresholds", m_TimingThresholds,
"Thresholds (MeV) to display ECL timing", std::vector<double> {5, 10, 50});
71 addParam(
"HitNumberUpperlimits", m_HitNumberUpperLimits,
72 "Upper limit (# of hits) to display hit multiplicity", std::vector<double> {10000, 1000, 700, 200});
73 addParam(
"WaveformOption", m_WaveformOption,
"Option (all,psd,logic,rand,dphy,other) to display waveform flow",
75 addParam(
"DPHYTTYP", m_DPHYTTYP,
76 "Flag to control trigger of delayed bhabha events; 0 - select events by 'bha_delay' trigger bit, 1 - select by TTYP_DPHY",
false);
86 TDirectory* oldDir = gDirectory;
90 TDirectory* dirDAQ =
dynamic_cast<TDirectory*
>(oldDir->Get(m_histogramDirectoryName.c_str()));
91 if (!dirDAQ) dirDAQ = oldDir->mkdir(m_histogramDirectoryName.c_str());
95 h_evtot =
new TH1F(
"event",
"Total event bank", 1, 0, 1);
96 h_evtot->SetOption(
"LIVE");
98 h_quality =
new TH1F(
"quality",
"Fit quality flag. 0-good, 1-integer overflow, 2-low amplitude, 3-bad chi2", 4, 0, 4);
99 h_quality->GetXaxis()->SetTitle(
"Flag number");
100 h_quality->SetFillColor(kPink - 4);
101 h_quality->SetOption(
"LIVE");
103 h_quality_other =
new TH1F(
"quality_other",
"Fit quality flag for unexpectedly saved waveforms", 4, 0, 4);
104 h_quality_other->GetXaxis()->SetTitle(
"Flag number. 0-good,1-int overflow,2-low amplitude,3-bad chi2");
105 h_quality_other->SetFillColor(kPink - 4);
106 h_quality_other->SetOption(
"LIVE");
108 h_bad_quality =
new TH1F(
"bad_quality",
"Fraction of hits with bad chi2 (qual=3) and E > 1 GeV vs Cell ID", 8736, 1, 8737);
109 h_bad_quality->GetXaxis()->SetTitle(
"Cell ID");
110 h_bad_quality->SetOption(
"LIVE");
112 h_trigtag1 =
new TH1F(
"trigtag1",
"Consistency b/w global event number and trigger tag. 0-good, 1-DQM error", 2, 0, 2);
113 h_trigtag1->GetXaxis()->SetTitle(
"Flag number");
114 h_trigtag1->SetOption(
"LIVE");
115 h_trigtag1->SetFillColor(kPink - 4);
117 h_adc_hits =
new TH1F(
"adc_hits",
"Fraction of high-energy hits (E > 50 MeV)", 1001, 0, 1.001);
118 h_adc_hits->GetXaxis()->SetTitle(
"Fraction");
119 h_adc_hits->SetOption(
"LIVE");
121 for (
const auto&
id : m_HitThresholds) {
122 std::string h_name, h_title;
123 h_name = str(boost::format(
"cid_Thr%1%MeV") %
id);
124 h_title = str(boost::format(
"Occupancy per Cell ID (E > %1% MeV)") %
id);
125 TH1F* h =
new TH1F(h_name.c_str(), h_title.c_str(), 8736, 1, 8737);
126 h->GetXaxis()->SetTitle(
"Cell ID");
127 h->GetYaxis()->SetTitle(
"Occupancy (hits / events_count)");
128 h->SetOption(
"LIVE");
132 for (
const auto&
id : m_TotalEnergyThresholds) {
133 std::string h_name, h_title;
134 h_name = str(boost::format(
"edep_Thr%1%MeV") %
id);
135 h_title = str(boost::format(
"Total energy (thr = %1% MeV)") %
id);
136 TH1F* h =
new TH1F(h_name.c_str(), h_title.c_str(), (
int)(100 * m_EnergyUpperThr), 0, m_EnergyUpperThr);
137 h->GetXaxis()->SetTitle(
"Energy, [GeV]");
138 h->SetOption(
"LIVE");
139 h_edeps.push_back(h);
142 for (
const auto&
id : m_TimingThresholds) {
143 std::string h_bar_name, h_bar_title;
144 std::string h_end_name, h_end_title;
145 h_bar_name = str(boost::format(
"time_barrel_Thr%1%MeV") %
id);
146 h_bar_title = str(boost::format(
"Reconstructed time for ECL barrel (E > %1% MeV)") %
id);
147 h_end_name = str(boost::format(
"time_endcaps_Thr%1%MeV") %
id);
148 h_end_title = str(boost::format(
"Reconstructed time for ECL endcaps (E > %1% MeV)") %
id);
149 TH1F* h_time_barrel =
new TH1F(h_bar_name.c_str(), h_bar_title.c_str(), 206, -1030, 1030);
150 TH1F* h_time_endcap =
new TH1F(h_end_name.c_str(), h_end_title.c_str(), 206, -1030, 1030);
151 h_time_barrel->GetXaxis()->SetTitle(
"Time, [ns]");
152 h_time_endcap->GetXaxis()->SetTitle(
"Time, [ns]");
153 h_time_barrel->SetOption(
"LIVE");
154 h_time_endcap->SetOption(
"LIVE");
155 h_time_barrels.push_back(h_time_barrel);
156 h_time_endcaps.push_back(h_time_endcap);
159 for (
const auto&
id : boost::combine(m_HitThresholds, m_HitNumberUpperLimits)) {
160 double id1 = 0, id2 = 0;
161 boost::tie(id1, id2) = id;
162 std::string h_name, h_title;
163 h_name = str(boost::format(
"ncev_Thr%1%MeV") % id1);
164 h_title = str(boost::format(
"Number of hits in event (E > %1% MeV)") % id1);
165 TH1F* h =
new TH1F(h_name.c_str(), h_title.c_str(), id2, 0, id2);
166 h->GetXaxis()->SetTitle(
"Number of hits");
167 h->SetOption(
"LIVE");
168 h_ncevs.push_back(h);
171 for (
int i = 0; i < ECL_CRATES; i++) {
173 std::string h_name, h_title;
174 h_name = str(boost::format(
"time_crate_%1%_Thr1GeV") % (crate));
175 h_title = str(boost::format(
"Reconstructed time for ECL crate #%1% with E > 1 GeV") % (crate));
176 TH1F* h =
new TH1F(h_name.c_str(), h_title.c_str(), 400, -100, 100);
177 h->GetXaxis()->SetTitle(
"Time [ns]");
178 h->SetOption(
"LIVE");
179 h_time_crate_Thr1GeV.push_back(h);
182 for (
const auto&
id : m_WaveformOption) {
183 if (
id !=
"all" &&
id !=
"psd" &&
id !=
"logic" &&
id !=
"rand" &&
id !=
"dphy" &&
id !=
"other")
184 B2WARNING(
"Waveform Options are not correctly assigned. They must be 'all', 'psd', 'logic', 'rand', 'dphy', 'other'!");
186 std::string h_cell_name;
187 if (
id ==
"other") h_title =
"Unexpectedly saved waveforms";
188 if (
id ==
"psd") h_title =
"#frac{Saved}{Expected} waveforms for high-energy hits (E > 50 MeV)";
189 if (
id ==
"logic") h_title =
"#frac{Saved}{Expected} waveforms for every 1000th event";
190 if (
id ==
"rand") h_title =
"#frac{Saved}{Expected} waveforms for random trigger events";
191 if (
id ==
"dphy") h_title =
"#frac{Saved}{Expected} waveforms for delayed bhabha (DPHY) events";
192 if (
id ==
"all") h_title =
"#frac{Saved}{Expected} waveforms for all events";
193 h_cell_name = str(boost::format(
"wf_cid_%1%") % (
id));
194 TH1F* h_cell =
new TH1F(h_cell_name.c_str(), h_title.c_str(), 8736, 1, 8737);
195 h_cell->GetXaxis()->SetTitle(
"Cell ID");
196 h_cell->SetOption(
"LIVE");
198 h_cell_psd_norm =
new TH1F(
"psd_cid",
"Normalization to psd hits for cid", 8736, 1, 8737);
199 h_cell_psd_norm->SetOption(
"LIVE");
202 h_evtot_logic =
new TH1F(
"event_logic",
"Event bank for logic", 1, 0, 1);
203 h_evtot_logic->SetOption(
"LIVE");
206 h_evtot_rand =
new TH1F(
"event_rand",
"Event bank for rand", 1, 0, 1);
207 h_evtot_rand->SetOption(
"LIVE");
210 h_evtot_dphy =
new TH1F(
"event_dphy",
"Event bank for dphy", 1, 0, 1);
211 h_evtot_dphy->SetOption(
"LIVE");
213 h_cells.push_back(h_cell);
218 h_trigtag2_trigid =
new TH2F(
"trigtag2_trigid",
"Internal data consistency vs crate. 0-good, 1-data corruption",
219 52, 1, 53, 11, -1, 10);
220 h_trigtag2_trigid->GetXaxis()->SetTitle(
"Crate ID (same as ECLCollector ID)");
221 h_trigtag2_trigid->GetYaxis()->SetTitle(
"Data consistency flag");
222 h_trigtag2_trigid->SetOption(
"LIVE");
224 h_pedmean_cellid =
new TProfile(
"pedmean_cellid",
"Pedestal vs Cell ID", 8736, 1, 8737);
225 h_pedmean_cellid->GetXaxis()->SetTitle(
"Cell ID");
226 h_pedmean_cellid->GetYaxis()->SetTitle(
"Ped. average (ADC units, #approx 0.05 MeV)");
227 h_pedmean_cellid->SetOption(
"LIVE");
229 h_pedrms_cellid =
new TProfile(
"pedrms_cellid",
"Pedestal stddev vs Cell ID",
231 h_pedrms_cellid->GetXaxis()->SetTitle(
"Cell ID");
232 h_pedrms_cellid->GetYaxis()->SetTitle(
"Ped. stddev (ADC units, #approx 0.05 MeV)");
233 h_pedrms_cellid->SetOption(
"LIVE");
235 h_pedrms_thetaid =
new TProfile(
"pedrms_thetaid",
"Pedestal stddev vs #theta ID",
237 h_pedrms_thetaid->GetXaxis()->SetTitle(
"#theta ID (0-12=FWD, 59-67=BWD endcap)");
238 h_pedrms_thetaid->GetYaxis()->SetTitle(
"Ped. stddev (ADC units, #approx 0.05 MeV)");
239 h_pedrms_thetaid->SetOption(
"LIVE");
241 h_trigtime_trigid =
new TH2F(
"trigtime_trigid",
"Trigger time vs Crate ID", 52, 1, 53, 145, 0, 145);
242 h_trigtime_trigid->GetXaxis()->SetTitle(
"Crate ID (same as ECLCollector ID)");
243 h_trigtime_trigid->GetYaxis()->SetTitle(
"Trigger time (only even, 0-142)");
244 h_trigtime_trigid->SetOption(
"LIVE");
254 m_ECLDigits.isRequired();
255 m_ECLCalDigits.isOptional();
256 m_ECLTrigs.isOptional();
257 m_ECLDsps.isOptional();
258 m_l1Trigger.isOptional();
260 if (!mapper.initFromDB()) B2FATAL(
"ECL DQM: Can't initialize eclChannelMapper");
262 ecltot.resize(m_TotalEnergyThresholds.size());
263 nhits.resize(m_HitNumberUpperLimits.size());
267 v_totalthrApsd.resize((m_calibrationThrApsd->getCalibVector()).size());
268 for (
size_t i = 0; i < v_totalthrApsd.size(); i++) v_totalthrApsd[i] = (
int)(m_calibrationThrApsd->getCalibVector())[i];
274 h_evtot_logic->Reset();
275 h_evtot_rand->Reset();
276 h_evtot_dphy->Reset();
278 h_quality_other->Reset();
279 h_bad_quality->Reset();
282 h_cell_psd_norm->Reset();
283 std::for_each(h_cids.begin(), h_cids.end(), [](
auto & it) {it->Reset();});
284 std::for_each(h_edeps.begin(), h_edeps.end(), [](
auto & it) {it->Reset();});
285 std::for_each(h_time_barrels.begin(), h_time_barrels.end(), [](
auto & it) {it->Reset();});
286 std::for_each(h_time_endcaps.begin(), h_time_endcaps.end(), [](
auto & it) {it->Reset();});
287 std::for_each(h_ncevs.begin(), h_ncevs.end(), [](
auto & it) {it->Reset();});
288 std::for_each(h_cells.begin(), h_cells.end(), [](
auto & it) {it->Reset();});
289 for (
int i = 0; i < ECL_CRATES; i++) h_time_crate_Thr1GeV[i]->Reset();
290 h_trigtag2_trigid->Reset();
291 h_pedmean_cellid->Reset();
292 h_pedrms_cellid->Reset();
293 h_pedrms_thetaid->Reset();
294 h_trigtime_trigid->Reset();
301 for (
auto& value : ecltot) value = 0;
302 for (
auto& value : nhits) value = 0;
303 bool bhatrig =
false;
305 if (m_l1Trigger.isValid() && m_DPHYTTYP) bhatrig = m_l1Trigger->getTimType() == TRGSummary::ETimingType::TTYP_DPHY;
306 else if (m_l1Trigger.isValid() && !m_DPHYTTYP) bhatrig = m_l1Trigger->testInput(
"bha_delay");
309 if (m_eventmetadata.isValid()) {
310 if (m_eventmetadata->getErrorFlag() != 0x10) {
311 m_iEvent = m_eventmetadata->getEvent();
313 for (
const auto&
id : m_WaveformOption) {
314 if (
id ==
"logic" && m_iEvent % 1000 == 999) h_evtot_logic->Fill(0);
315 if (
id ==
"rand" && m_l1Trigger.isValid() &&
316 m_l1Trigger->getTimType() == TRGSummary::ETimingType::TTYP_RAND) h_evtot_rand->Fill(0);
317 if (
id ==
"dphy" && m_l1Trigger.isValid() &&
318 bhatrig) h_evtot_dphy->Fill(0);
323 for (
auto& aECLDigit : m_ECLDigits) {
324 int i = aECLDigit.getCellId() - 1;
325 h_quality->Fill(aECLDigit.getQuality());
326 if (aECLDigit.getAmp() > 2.e04 && aECLDigit.getQuality() == 3) h_bad_quality->Fill(aECLDigit.getCellId());
327 if (aECLDigit.getAmp() >= (v_totalthrApsd[i] / 4 * 4)) NDigits ++;
328 for (
const auto&
id : m_WaveformOption) {
329 if (
id !=
"psd")
continue;
330 else if (
id ==
"psd" && (m_iEvent % 1000 == 999 ||
331 (m_l1Trigger.isValid() && m_l1Trigger->getTimType() == TRGSummary::ETimingType::TTYP_RAND) ||
332 (m_l1Trigger.isValid() && bhatrig) ||
333 aECLDigit.getAmp() < (v_totalthrApsd[i] / 4 * 4)))
continue;
334 h_cell_psd_norm->Fill(aECLDigit.getCellId());
338 for (
auto& aECLTrig : m_ECLTrigs) {
339 double itrg = aECLTrig.getTimeTrig();
342 int tg = (int)itrg - 2 * ((
int)itrg / 8);
343 h_trigtime_trigid->Fill(aECLTrig.getTrigId(), tg);
344 trigtag1 += aECLTrig.getTrigTag();
345 h_trigtag2_trigid->Fill(aECLTrig.getTrigId(), aECLTrig.getTrigTagQualityFlag());
348 if (m_ECLTrigs.getEntries() > 0) {
350 trigtag1 /= m_ECLTrigs.getEntries();
351 int compar = (65535 & m_iEvent);
352 if (compar == trigtag1) flagtag = 0;
353 h_trigtag1->Fill(flagtag);
356 for (
auto& aECLCalDigit : m_ECLCalDigits) {
357 int cid = aECLCalDigit.getCellId();
358 double energy = aECLCalDigit.getEnergy();
359 double timing = aECLCalDigit.getTime();
361 for (
const auto&
id : m_HitThresholds) {
362 auto scale =
id / 1000.;
363 auto index = std::distance(m_HitThresholds.begin(), std::find(m_HitThresholds.begin(), m_HitThresholds.end(),
id));
364 if (energy > scale) {
365 h_cids[index]->Fill(cid);
370 for (
const auto&
id : m_TotalEnergyThresholds) {
371 auto scale =
id / 1000.;
372 auto index = std::distance(m_TotalEnergyThresholds.begin(), std::find(m_TotalEnergyThresholds.begin(),
373 m_TotalEnergyThresholds.end(),
id));
374 if (energy > scale) ecltot[index] += energy;
377 for (
const auto&
id : m_TimingThresholds) {
378 auto scale =
id / 1000.;
379 auto index = std::distance(m_TimingThresholds.begin(), std::find(m_TimingThresholds.begin(), m_TimingThresholds.end(),
id));
380 if (energy > scale) {
381 if (cid > ECL_FWD_CHANNELS && cid <= ECL_FWD_CHANNELS + ECL_BARREL_CHANNELS) h_time_barrels[index]->Fill(timing);
382 else h_time_endcaps[index]->Fill(timing);
386 if (energy > 1.000) h_time_crate_Thr1GeV[mapper.getCrateID(cid) - 1]->Fill(timing);
389 for (
auto& h : h_edeps) {
390 auto index = std::distance(h_edeps.begin(), std::find(h_edeps.begin(), h_edeps.end(), h));
391 h->Fill(ecltot[index]);
394 for (
auto& h : h_ncevs) {
395 auto index = std::distance(h_ncevs.begin(), std::find(h_ncevs.begin(), h_ncevs.end(), h));
396 h->Fill(nhits[index]);
399 for (
auto& aECLDsp : m_ECLDsps) {
400 int i = aECLDsp.getCellId() - 1;
401 aECLDsp.getDspA(m_DspArray[i]);
402 m_PedestalMean[i] = 0;
403 m_PedestalRms[i] = 0;
405 for (
int j = 0; j < 16; j++) m_PedestalMean[i] += m_DspArray[i][j];
406 m_PedestalMean[i] /= 16;
407 h_pedmean_cellid->Fill(aECLDsp.getCellId(), m_PedestalMean[i]);
409 for (
int j = 0; j < 16; j++) m_PedestalRms[i] += pow(m_DspArray[i][j] - m_PedestalMean[i], 2);
410 m_PedestalRms[i] = sqrt(m_PedestalRms[i] / 15.);
411 h_pedrms_cellid->Fill(aECLDsp.getCellId(), m_PedestalRms[i]);
413 h_pedrms_thetaid->Fill(m_geom->GetThetaID(), m_PedestalRms[i]);
417 for (
const auto&
id : m_WaveformOption) {
418 auto index = std::distance(m_WaveformOption.begin(), std::find(m_WaveformOption.begin(), m_WaveformOption.end(),
id));
419 if (
id !=
"all" &&
id !=
"psd" &&
id !=
"logic" &&
id !=
"rand" &&
id !=
"dphy" &&
id !=
"other")
continue;
420 else if (
id ==
"psd" && (m_iEvent % 1000 == 999 ||
421 (m_l1Trigger.isValid() && m_l1Trigger->getTimType() == TRGSummary::ETimingType::TTYP_RAND) ||
422 (m_l1Trigger.isValid() && bhatrig) ||
423 !aECLDigit || aECLDigit->
getAmp() < (v_totalthrApsd[i] / 4 * 4)))
continue;
424 else if (
id ==
"logic" && m_iEvent % 1000 != 999)
continue;
425 else if (
id ==
"rand" && (m_iEvent % 1000 == 999 || !m_l1Trigger.isValid() ||
426 m_l1Trigger->getTimType() != TRGSummary::ETimingType::TTYP_RAND))
continue;
427 else if (
id ==
"dphy" && (m_iEvent % 1000 == 999 || !m_l1Trigger.isValid() ||
429 else if (
id ==
"other" && (m_iEvent % 1000 == 999 ||
430 (m_l1Trigger.isValid() && m_l1Trigger->getTimType() == TRGSummary::ETimingType::TTYP_RAND) ||
431 (m_l1Trigger.isValid() && bhatrig) ||
432 (aECLDigit && aECLDigit->
getAmp() >= (v_totalthrApsd[i] / 4 * 4))))
continue;
433 h_cells[index]->Fill(aECLDsp.getCellId());
434 if (
id ==
"other" && aECLDigit) h_quality_other->Fill(aECLDigit->
getQuality());
437 if (m_ECLDigits.getEntries() > 0)
438 h_adc_hits->Fill((
double)NDigits / (
double)m_ECLDigits.getEntries());