12 #include <framework/core/HistoModule.h>
15 #include <top/modules/TOPDQM/TOPDQMModule.h>
16 #include <top/geometry/TOPGeometryPar.h>
19 #include <framework/datastore/StoreArray.h>
20 #include <framework/datastore/StoreObjPtr.h>
23 #include <framework/gearbox/Unit.h>
24 #include <framework/gearbox/Const.h>
25 #include <framework/logging/Logger.h>
28 #include <top/dataobjects/TOPDigit.h>
29 #include <top/dataobjects/TOPPull.h>
30 #include <top/dataobjects/TOPLikelihood.h>
31 #include <mdst/dataobjects/Track.h>
35 #include "TDirectory.h"
38 #include <boost/format.hpp>
64 setDescription(
"TOP DQM histogrammer");
65 setPropertyFlags(c_ParallelProcessingCertified);
68 addParam(
"histogramDirectoryName", m_histogramDirectoryName,
69 "histogram directory in ROOT file",
string(
"TOP"));
70 addParam(
"momentumCut", m_momentumCut,
71 "momentum cut used to histogram pulls etc.", 2.0);
72 addParam(
"pValueCut", m_pValueCut,
73 "track p-value cut used to histogram pulls etc.", 0.001);
74 addParam(
"usePionID", m_usePionID,
75 "use pion ID from TOP to histogram pulls etc.",
true);
76 addParam(
"cutNphot", m_cutNphot,
"Cut on total number of photons", 3);
80 TOPDQMModule::~TOPDQMModule()
84 void TOPDQMModule::defineHisto()
87 TDirectory* oldDir = gDirectory;
88 oldDir->mkdir(m_histogramDirectoryName.c_str())->cd();
91 const auto* geo = TOPGeometryPar::Instance()->getGeometry();
92 m_numModules = geo->getNumModules();
93 double bunchTimeSep = geo->getNominalTDC().getSyncTimeBase() / 24;
95 m_BoolEvtMonitor =
new TH1F(
"BoolEvtMonitor",
"Event desynchronization monitoring",
97 m_BoolEvtMonitor->GetXaxis()->SetTitle(
"good/bad event entries");
99 m_recoTime =
new TH1F(
"recoTime",
"reco: time distribution",
101 m_recoTime->GetXaxis()->SetTitle(
"time [ns]");
103 m_recoTimeBg =
new TH1F(
"recoTimeBg",
"reco: time distribution (bkg)",
105 m_recoTimeBg->GetXaxis()->SetTitle(
"time [ns]");
107 m_recoTimeMinT0 =
new TH1F(
"recoTimeMinT0",
"reco: time in respect to first PDF peak",
109 m_recoTimeMinT0->GetXaxis()->SetTitle(
"time [ns]");
111 m_recoTimeDiff =
new TH1F(
"recoTimeDiff",
"reco: time resolution",
113 m_recoTimeDiff->GetXaxis()->SetTitle(
"time residual [ns]");
115 m_recoPull =
new TH1F(
"recoPull",
"reco: pulls",
117 m_recoPull->GetXaxis()->SetTitle(
"pull");
119 m_recoTimeDiff_Phic =
new TH2F(
"recoTimeDiff_Phic",
120 "reco: time resolution vs. phiCer",
121 90, -180, 180, 100, -1.0, 1.0);
122 m_recoTimeDiff_Phic->GetXaxis()->SetTitle(
"Cerenkov azimuthal angle [deg]");
123 m_recoTimeDiff_Phic->GetYaxis()->SetTitle(
"time residuals [ns]");
125 m_recoPull_Phic =
new TProfile(
"recoPull_Phic",
126 "reco: pulls vs phiCer",
127 90, -180, 180, -10, 10,
"S");
128 m_recoPull_Phic->GetXaxis()->SetTitle(
"Cerenkov azimuthal angle [deg]");
129 m_recoPull_Phic->GetYaxis()->SetTitle(
"pulls");
132 m_goodHits =
new TH1F(
"goodHits",
"Number of good hits per bar", m_numModules, 0.5, m_numModules + 0.5);
133 m_badHits =
new TH1F(
"badHits",
"Number of bad hits per bar", m_numModules, 0.5, m_numModules + 0.5);
134 m_goodHits->SetOption(
"LIVE");
135 m_badHits->SetOption(
"LIVE");
136 m_goodHits->SetMinimum(0);
137 m_badHits->SetMinimum(0);
139 m_goodHits->GetXaxis()->SetTitle(
"slot no.");
140 m_goodHits->GetYaxis()->SetTitle(
"hits / slot");
141 m_badHits->GetXaxis()->SetTitle(
"slot number");
142 m_badHits->GetYaxis()->SetTitle(
"hits / slot");
145 m_window_vs_slot =
new TH2F(
"window_vs_slot",
"Distribution of hits: raw timing",
146 16, 0.5, 16.5, 512, 0, 512);
147 m_window_vs_slot->SetXTitle(
"slot number");
148 m_window_vs_slot->SetYTitle(
"window number w.r.t reference window");
149 m_window_vs_slot->SetOption(
"LIVE");
150 m_window_vs_slot->SetStats(kFALSE);
152 m_bunchOffset =
new TH1F(
"bunchOffset",
"Reconstructed bunch: current offset",
153 100, -bunchTimeSep / 2, bunchTimeSep / 2);
154 m_bunchOffset->SetXTitle(
"offset [ns]");
155 m_bunchOffset->SetYTitle(
"events/bin");
156 m_bunchOffset->SetOption(
"LIVE");
157 m_bunchOffset->SetMinimum(0);
159 m_time =
new TH1F(
"goodHitTimes",
"Time distribution of good hits",
161 m_time->SetXTitle(
"time [ns]");
162 m_time->SetYTitle(
"hits/bin");
163 m_time->SetOption(
"LIVE");
164 m_time->SetMinimum(0);
167 m_goodHitsPerEventAll =
new TH1F(
"goodHitsPerEventAll",
"Number of good hits per event", MaxEvents, 0, MaxEvents);
168 m_badHitsPerEventAll =
new TH1F(
"badHitsPerEventAll",
"Number of bad hits per event", MaxEvents, 0, MaxEvents);
169 m_goodHitsPerEventAll->SetOption(
"LIVE");
170 m_badHitsPerEventAll->SetOption(
"LIVE");
171 m_goodHitsPerEventAll->SetMinimum(0);
172 m_badHitsPerEventAll->SetMinimum(0);
173 m_goodHitsPerEventAll->GetXaxis()->SetTitle(
"hits / event");
174 m_goodHitsPerEventAll->GetYaxis()->SetTitle(
"Events");
175 m_badHitsPerEventAll->GetXaxis()->SetTitle(
"hits / event");
176 m_badHitsPerEventAll->GetYaxis()->SetTitle(
"Events");
180 m_goodTDCAll =
new TH1F(
"goodTDCAll",
"Raw time distribution of good hits",
181 BinNumRT, 0, MaxRawTime);
182 m_goodTDCAll->SetXTitle(
"raw time [samples]");
183 m_goodTDCAll->SetYTitle(
"hits / sample");
184 m_goodTDCAll->SetOption(
"LIVE");
185 m_goodTDCAll->SetMinimum(0);
187 m_badTDCAll =
new TH1F(
"badTDCAll",
"Raw time distribution of bad hits",
188 BinNumRT, 0, MaxRawTime);
189 m_badTDCAll->SetXTitle(
"raw time [samples]");
190 m_badTDCAll->SetYTitle(
"hits / sample");
191 m_badTDCAll->SetOption(
"LIVE");
192 m_badTDCAll->SetMinimum(0);
194 m_goodHitsPerEventProf =
new TProfile(
"goodHitsPerEventProf",
"Good hits per event vs. slot number",
195 16, 0.5, 16.5, 0, MaxEvents);
196 m_goodHitsPerEventProf->SetXTitle(
"slot number");
197 m_goodHitsPerEventProf->SetYTitle(
"hits per event");
198 m_goodHitsPerEventProf->SetOption(
"LIVE");
199 m_goodHitsPerEventProf->SetStats(kFALSE);
200 m_goodHitsPerEventProf->SetMinimum(0);
202 m_badHitsPerEventProf =
new TProfile(
"badHitsPerEventProf",
"Bad hits per event vs. slot number",
203 16, 0.5, 16.5, 0, MaxEvents);
204 m_badHitsPerEventProf->SetXTitle(
"slot number");
205 m_badHitsPerEventProf->SetYTitle(
"hits per event");
206 m_badHitsPerEventProf->SetOption(
"LIVE");
207 m_badHitsPerEventProf->SetStats(kFALSE);
208 m_badHitsPerEventProf->SetMinimum(0);
210 m_TOPOccAfterInjLER =
new TH1F(
"TOPOccInjLER",
"TOPOccInjLER/Time;Time in #mus;Nhits/Time (#mus bins)", 4000, 0, 20000);
211 m_TOPOccAfterInjHER =
new TH1F(
"TOPOccInjHER",
"TOPOccInjHER/Time;Time in #mus;Nhits/Time (#mus bins)", 4000, 0, 20000);
212 m_TOPEOccAfterInjLER =
new TH1F(
"TOPEOccInjLER",
"TOPEOccInjLER/Time;Time in #mus;Triggers/Time (#mus bins)", 4000, 0, 20000);
213 m_TOPEOccAfterInjHER =
new TH1F(
"TOPEOccInjHER",
"TOPEOccInjHER/Time;Time in #mus;Triggers/Time (#mus bins)", 4000, 0, 20000);
215 for (
int i = 0; i < m_numModules; i++) {
222 name = str(format(
"window_vs_asic_%1%") % (module));
223 title = str(format(
"Distribution of hits: raw timing for slot #%1%") % (module));
224 h2 =
new TH2F(name.c_str(), title.c_str(), 64, 0, 64, 512, 0, 512);
225 h2->SetOption(
"LIVE");
226 h2->SetStats(kFALSE);
227 h2->SetXTitle(
"ASIC number");
228 h2->SetYTitle(
"window number w.r.t reference window");
230 m_window_vs_asic.push_back(h2);
232 name = str(format(
"good_hits_xy_%1%") % (module));
233 title = str(format(
"Number of good hits in x-y for slot #%1%") % (module));
234 h2 =
new TH2F(name.c_str(), title.c_str(), 64, 0.5, 64.5, 8, 0.5, 8.5);
235 h2->SetOption(
"LIVE");
236 h2->SetStats(kFALSE);
237 h2->GetXaxis()->SetTitle(
"pixel column");
238 h2->GetYaxis()->SetTitle(
"pixel row");
240 m_goodHitsXY.push_back(h2);
242 name = str(format(
"bad_hits_xy_%1%") % (module));
243 title = str(format(
"Number of bad hits in x-y for slot #%1%") % (module));
244 h2 =
new TH2F(name.c_str(), title.c_str(), 64, 0.5, 64.5, 8, 0.5, 8.5);
245 h2->SetOption(
"LIVE");
246 h2->SetStats(kFALSE);
247 h2->GetXaxis()->SetTitle(
"pixel column");
248 h2->GetYaxis()->SetTitle(
"pixel row");
250 m_badHitsXY.push_back(h2);
252 name = str(format(
"good_hits_asics_%1%") % (module));
253 title = str(format(
"Number of good hits for asics for slot #%1%") % (module));
254 h2 =
new TH2F(name.c_str(), title.c_str(), 64, 0, 64, 8, 0, 8);
255 h2->SetOption(
"LIVE");
256 h2->SetStats(kFALSE);
257 h2->GetXaxis()->SetTitle(
"ASIC number");
258 h2->GetYaxis()->SetTitle(
"ASIC channel");
260 m_goodHitsAsics.push_back(h2);
262 name = str(format(
"bad_hits_asics_%1%") % (module));
263 title = str(format(
"Number of bad hits for asics for slot #%1%") % (module));
264 h2 =
new TH2F(name.c_str(), title.c_str(), 64, 0, 64, 8, 0, 8);
265 h2->SetOption(
"LIVE");
266 h2->SetStats(kFALSE);
267 h2->GetXaxis()->SetTitle(
"ASIC number");
268 h2->GetYaxis()->SetTitle(
"ASIC channel");
270 m_badHitsAsics.push_back(h2);
272 name = str(format(
"good_TDC_%1%") % (module));
273 title = str(format(
"Raw time distribution of good hits for slot #%1%") % (module));
274 h1 =
new TH1F(name.c_str(), title.c_str(), BinNumRT, 0, MaxRawTime);
275 h1->SetOption(
"LIVE");
276 h1->GetXaxis()->SetTitle(
"raw time [samples]");
277 h1->GetYaxis()->SetTitle(
"hits per sample");
279 m_goodTdc.push_back(h1);
281 name = str(format(
"bad_TDC_%1%") % (module));
282 title = str(format(
"Raw time distribution of bad hits for slot #%1%") % (module));
283 h1 =
new TH1F(name.c_str(), title.c_str(), BinNumRT, 0, MaxRawTime);
284 h1->SetOption(
"LIVE");
285 h1->GetXaxis()->SetTitle(
"raw time [samples]");
286 h1->GetYaxis()->SetTitle(
"hits per sample");
288 m_badTdc.push_back(h1);
290 name = str(format(
"good_timing_%1%") % (module));
291 title = str(format(
"Timing distribution of good hits for slot #%1%") % (module));
292 h1 =
new TH1F(name.c_str(), title.c_str(), 100, -20, 80);
293 h1->SetOption(
"LIVE");
294 h1->GetXaxis()->SetTitle(
"time [ns]");
295 h1->GetYaxis()->SetTitle(
"hits per time bin");
297 m_goodTiming.push_back(h1);
299 name = str(format(
"good_channel_hits_%1%") % (module));
300 title = str(format(
"Number of good hits by channel for slot #%1%") % (module));
301 int numPixels = geo->getModule(i + 1).getPMTArray().getNumPixels();
302 h1 =
new TH1F(name.c_str(), title.c_str(), numPixels, 0, numPixels);
303 h1->SetOption(
"LIVE");
304 h1->GetXaxis()->SetTitle(
"channel number");
305 h1->GetYaxis()->SetTitle(
"hits pre channel");
307 m_goodChannelHits.push_back(h1);
309 name = str(format(
"bad_channel_hits_%1%") % (module));
310 title = str(format(
"Number of bad hits by channel for slot #%1%") % (module));
311 h1 =
new TH1F(name.c_str(), title.c_str(), numPixels, 0, numPixels);
312 h1->SetOption(
"LIVE");
313 h1->GetXaxis()->SetTitle(
"channel number");
314 h1->GetYaxis()->SetTitle(
"hits pre channel");
316 m_badChannelHits.push_back(h1);
318 name = str(format(
"good_hits_per_event%1%") % (module));
319 title = str(format(
"Number of good hits per event for slot #%1%") % (module));
320 h1 =
new TH1F(name.c_str(), title.c_str(), MaxEvents, 0, MaxEvents);
321 h1->SetOption(
"LIVE");
322 h1->GetXaxis()->SetTitle(
"hits / event");
324 m_goodHitsPerEvent.push_back(h1);
326 name = str(format(
"bad_hits_per_event%1%") % (module));
327 title = str(format(
"Number of bad hits per event for slot #%1%") % (module));
328 h1 =
new TH1F(name.c_str(), title.c_str(), MaxEvents, 0, MaxEvents);
329 h1->SetOption(
"LIVE");
330 h1->GetXaxis()->SetTitle(
"hits / event");
332 m_badHitsPerEvent.push_back(h1);
334 name = str(format(
"good_hits_xy_track_%1%") % (module));
335 title = str(format(
"Hits per track, each channel, slot #%1%") % (module));
336 h3 =
new TProfile2D(name.c_str(), title.c_str(), 64, 0.5, 64.5, 8, 0.5, 8.5, 0, 1000);
337 h3->SetOption(
"LIVE");
338 h3->SetStats(kFALSE);
339 h3->GetXaxis()->SetTitle(
"pixel column");
340 h3->GetYaxis()->SetTitle(
"pixel row");
342 m_goodHitsXYTrack.push_back(h3);
344 name = str(format(
"good_hits_xy_track_bkg_%1%") % (module));
345 title = str(format(
"Hits per bkg track, each channel, slot #%1%") % (module));
346 h3 =
new TProfile2D(name.c_str(), title.c_str(), 64, 0.5, 64.5, 8, 0.5, 8.5, 0, 1000);
347 h3->SetOption(
"LIVE");
348 h3->SetStats(kFALSE);
349 h3->GetXaxis()->SetTitle(
"pixel column");
350 h3->GetYaxis()->SetTitle(
"pixel row");
352 m_goodHitsXYTrackBkg.push_back(h3);
359 void TOPDQMModule::initialize()
365 m_rawFTSW.isOptional();
366 m_digits.isRequired();
367 m_recBunch.isOptional();
368 m_tracks.isOptional();
372 void TOPDQMModule::beginRun()
374 m_BoolEvtMonitor->Reset();
376 m_recoTimeDiff->Reset();
377 m_recoTimeDiff_Phic->Reset();
379 m_recoPull_Phic->Reset();
381 m_recoTimeBg->Reset();
382 m_recoTimeMinT0->Reset();
386 m_window_vs_slot->Reset();
387 m_bunchOffset->Reset();
389 m_goodTDCAll->Reset();
390 m_badTDCAll->Reset();
391 m_goodHitsPerEventProf->Reset();
392 m_goodHitsPerEventAll->Reset();
393 m_badHitsPerEventProf->Reset();
394 m_badHitsPerEventAll->Reset();
395 m_TOPOccAfterInjLER->Reset();
396 m_TOPOccAfterInjHER->Reset();
397 m_TOPEOccAfterInjLER->Reset();
398 m_TOPEOccAfterInjHER->Reset();
400 for (
int i = 0; i < m_numModules; i++) {
401 m_window_vs_asic[i]->Reset();
402 m_goodHitsXY[i]->Reset();
403 m_badHitsXY[i]->Reset();
404 m_goodHitsAsics[i]->Reset();
405 m_badHitsAsics[i]->Reset();
406 m_goodTdc[i]->Reset();
407 m_badTdc[i]->Reset();
408 m_goodTiming[i]->Reset();
409 m_goodChannelHits[i]->Reset();
410 m_badChannelHits[i]->Reset();
411 m_goodHitsPerEvent[i]->Reset();
412 m_badHitsPerEvent[i]->Reset();
413 m_goodHitsXYTrack[i]->Reset();
414 m_goodHitsXYTrackBkg[i]->Reset();
418 void TOPDQMModule::event()
421 bool recBunchValid =
false;
422 if (m_recBunch.isValid()) {
423 recBunchValid = m_recBunch->isReconstructed();
427 m_bunchOffset->Fill(m_recBunch->getCurrentOffset());
430 std::vector<int> n_good(16, 0);
431 std::vector<int> n_bad(16, 0);
432 std::vector<int> n_good_first(16, 0);
433 std::vector<int> n_good_second(16, 0);
434 std::vector<int> n_good_pixel_hits(16 * 512, 0);
436 int Ndigits = m_digits.getEntries();
438 for (
const auto& digit : m_digits) {
439 int x = digit.getFirstWindow() != m_digits[0]->getFirstWindow() ? 1 : 0 ;
440 m_BoolEvtMonitor->Fill(x);
444 for (
const auto& digit : m_digits) {
445 int i = digit.getModuleID() - 1;
446 if (i < 0 || i >= m_numModules) {
447 B2ERROR(
"Invalid module ID found in TOPDigits: ID = " << i + 1);
451 int ch = digit.getChannel();
452 int asic_no = ch / 8, asic_ch = ch % 8;
454 m_window_vs_slot->Fill(digit.getModuleID(), digit.getRawTime() / 64 + 220);
455 m_window_vs_asic[i]->Fill(digit.getChannel() / 8, digit.getRawTime() / 64 + 220);
457 if (digit.getHitQuality() != TOPDigit::c_Junk) {
458 m_goodHits->Fill(i + 1);
459 m_goodHitsXY[i]->Fill(digit.getPixelCol(), digit.getPixelRow());
460 m_goodHitsAsics[i]->Fill(asic_no, asic_ch);
461 m_goodTdc[i]->Fill(digit.getRawTime());
462 m_goodTDCAll->Fill(digit.getRawTime());
464 m_goodTiming[i]->Fill(digit.getTime());
465 m_time->Fill(digit.getTime());
467 m_goodChannelHits[i]->Fill(digit.getChannel());
468 if (digit.getTime() > 0 && digit.getTime() < 20) n_good_first[i]++;
469 if (digit.getTime() > 20 && digit.getTime() < 50) n_good_second[i]++;
470 n_good_pixel_hits[(digit.getModuleID() - 1) * 512 + (digit.getPixelID() - 1)]++;
473 m_badHits->Fill(i + 1);
474 m_badHitsXY[i]->Fill(digit.getPixelCol(), digit.getPixelRow());
475 m_badHitsAsics[i]->Fill(asic_no, asic_ch);
476 m_badTdc[i]->Fill(digit.getRawTime());
477 m_badTDCAll->Fill(digit.getRawTime());
478 m_badChannelHits[i]->Fill(digit.getChannel());
483 for (
int i = 0; i < 16; i++) {
484 m_goodHitsPerEventProf->Fill(i + 1, n_good[i]);
485 m_badHitsPerEventProf->Fill(i + 1, n_bad[i]);
486 m_goodHitsPerEvent[i]->Fill(n_good[i]);
487 m_goodHitsPerEventAll->Fill(n_good[i]);
488 m_badHitsPerEvent[i]->Fill(n_bad[i]);
489 m_badHitsPerEventAll->Fill(n_bad[i]);
491 bool slot_has_track = (n_good_first[i] + n_good_second[i]) > m_cutNphot;
492 for (
int j = 0; j < 512; j++) {
493 int col = j % 64 + 1, row = j / 64 + 1;
495 m_goodHitsXYTrack[i]->Fill(col, row, n_good_pixel_hits[i * 512 + j]);
497 m_goodHitsXYTrackBkg[i]->Fill(col, row, n_good_pixel_hits[i * 512 + j]);
501 for (
const auto& track : m_tracks) {
502 const auto* trackFit = track.getTrackFitResultWithClosestMass(Const::pion);
503 if (!trackFit)
continue;
504 if (trackFit->getMomentum().Mag() < m_momentumCut)
continue;
505 if (trackFit->getPValue() < m_pValueCut)
continue;
509 if (top->getLogL_pi() < top->getLogL_K())
continue;
510 if (top->getLogL_pi() < top->getLogL_p())
continue;
513 const auto pulls = track.getRelationsWith<
TOPPull>();
514 for (
const auto& pull : pulls) {
515 if (pull.isSignal()) {
516 double phiCer = pull.
getPhiCer() / Unit::deg;
517 m_recoTimeDiff->Fill(pull.getTimeDiff(), pull.getWeight());
518 m_recoTimeDiff_Phic->Fill(phiCer, pull.getTimeDiff(), pull.getWeight());
519 m_recoPull->Fill(pull.getPull(), pull.getWeight());
520 m_recoPull_Phic->Fill(phiCer, pull.getPull(), pull.getWeight());
522 m_recoTime->Fill(pull.getTime());
523 m_recoTimeBg->Fill(pull.getTime(), pull.getWeight());
524 m_recoTimeMinT0->Fill(pull.getTimeDiff());
529 for (
auto& it : m_rawFTSW) {
530 B2DEBUG(29,
"TTD FTSW : " << hex << it.GetTTUtime(0) <<
" " << it.GetTTCtime(0) <<
" EvtNr " << it.GetEveNo(0) <<
" Type " <<
531 (it.GetTTCtimeTRGType(0) & 0xF) <<
" TimeSincePrev " << it.GetTimeSincePrevTrigger(0) <<
" TimeSinceInj " <<
532 it.GetTimeSinceLastInjection(0) <<
" IsHER " << it.GetIsHER(0) <<
" Bunch " << it.GetBunchNumber(0));
533 auto difference = it.GetTimeSinceLastInjection(0);
534 if (difference != 0x7FFFFFFF) {
535 unsigned int nentries = m_digits.getEntries();
536 float diff2 = difference / 127.;
537 if (it.GetIsHER(0)) {
538 m_TOPOccAfterInjHER->Fill(diff2, nentries);
539 m_TOPEOccAfterInjHER->Fill(diff2);
541 m_TOPOccAfterInjLER->Fill(diff2, nentries);
542 m_TOPEOccAfterInjLER->Fill(diff2);
550 void TOPDQMModule::endRun()
554 void TOPDQMModule::terminate()