10#include <ecl/modules/eclDQM/eclDQM.h>
13#include <ecl/dataobjects/ECLCalDigit.h>
14#include <ecl/dataobjects/ECLDigit.h>
15#include <ecl/dataobjects/ECLDsp.h>
16#include <ecl/dataobjects/ECLElementNumbers.h>
17#include <ecl/dataobjects/ECLTrig.h>
18#include <ecl/geometry/ECLGeometryPar.h>
19#include <ecl/mapper/ECLChannelMapper.h>
22#include <analysis/dataobjects/ParticleList.h>
23#include <analysis/VariableManager/Manager.h>
24#include <analysis/variables/ECLVariables.h>
25#include <framework/core/HistoModule.h>
26#include <framework/dataobjects/EventMetaData.h>
27#include <framework/gearbox/Unit.h>
28#include <framework/logging/Logger.h>
29#include <mdst/dataobjects/TRGSummary.h>
32#include <boost/format.hpp>
33#include <boost/range/combine.hpp>
34#include <boost/range/adaptor/indexed.hpp>
37#include <TDirectory.h>
52using namespace boost::adaptors;
62 "Primary module for ECL Data Quality Monitor.\n"
63 "This module provides a large set of low-level histograms: occupancy, time distribution, number of saved waveforms, etc.");
69 "histogram directory in ROOT file", std::string(
"ECL"));
71 addParam(
"HitThresholds",
m_HitThresholds,
"Thresholds to display hit occupancy, MeV", std::vector<double> {0, 5, 10, 50});
75 "Upper limit (# of hits) to display hit multiplicity", std::vector<double> {10000, 1000, 700, 200});
79 "Flag to control trigger of delayed bhabha events; 0 - select events by 'bha_delay' trigger bit, 1 - select by TTYP_DPHY",
false);
90 TDirectory* oldDir = gDirectory;
99 h_evtot =
new TH1F(
"event",
"Total event bank", 1, 0, 1);
101 h_quality =
new TH1F(
"quality",
"Fit quality flag. 0-good, 1-integer overflow, 2-low amplitude, 3-bad chi2", 4, 0, 4);
102 h_quality->GetXaxis()->SetTitle(
"Quality flag");
103 h_quality->GetYaxis()->SetTitle(
"ECL hits count");
106 h_quality_other =
new TH1F(
"quality_other",
"Fit quality flag for unexpectedly saved waveforms", 4, 0, 4);
107 h_quality_other->GetXaxis()->SetTitle(
"Quality flag. 0-good,1-int overflow,2-low amplitude,3-bad chi2");
110 h_bad_quality =
new TH1F(
"bad_quality",
"Fraction of hits with bad chi2 (qual=3) and E > 1 GeV vs Cell ID",
115 h_trigtag1 =
new TH1F(
"trigtag1",
"Consistency b/w global event number and trigger tag. 0-good, 1-DQM error", 2, 0, 2);
116 h_trigtag1->GetXaxis()->SetTitle(
"Flag value");
117 h_trigtag1->GetYaxis()->SetTitle(
"Events count");
121 h_adc_hits =
new TH1F(
"adc_hits",
"Fraction of high-energy hits (E > 50 MeV)", 1001, 0, 1.001);
123 h_adc_hits->GetYaxis()->SetTitle(
"Events count");
126 "Number of hits with timing outside #pm 100 ns per Crate ID (E > 1 GeV)",
132 std::string h_name, h_title;
133 h_name = str(boost::format(
"cid_Thr%1%MeV") %
id);
134 h_title = str(boost::format(
"Occupancy per Cell ID (E > %1% MeV)") %
id);
135 TH1F* h =
new TH1F(h_name.c_str(), h_title.c_str(),
137 h->GetXaxis()->SetTitle(
"Cell ID");
138 h->GetYaxis()->SetTitle(
"Occupancy (hits / events_count)");
143 std::string h_name, h_title;
144 h_name = str(boost::format(
"edep_Thr%1%MeV") %
id);
145 h_title = str(boost::format(
"Total energy (thr = %1% MeV)") %
id);
147 h->GetXaxis()->SetTitle(
"Energy, [GeV]");
152 std::string h_bar_name, h_bar_title;
153 std::string h_end_name, h_end_title;
154 h_bar_name = str(boost::format(
"time_barrel_Thr%1%MeV") %
id);
155 h_bar_title = str(boost::format(
"Reconstructed time for ECL barrel (E > %1% MeV)") %
id);
156 h_end_name = str(boost::format(
"time_endcaps_Thr%1%MeV") %
id);
157 h_end_title = str(boost::format(
"Reconstructed time for ECL endcaps (E > %1% MeV)") %
id);
158 TH1F* h_time_barrel =
new TH1F(h_bar_name.c_str(), h_bar_title.c_str(), 206, -1030, 1030);
159 TH1F* h_time_endcap =
new TH1F(h_end_name.c_str(), h_end_title.c_str(), 206, -1030, 1030);
160 h_time_barrel->GetXaxis()->SetTitle(
"Time, [ns]");
161 h_time_endcap->GetXaxis()->SetTitle(
"Time, [ns]");
167 double id1 = 0, id2 = 0;
168 boost::tie(id1, id2) = id;
169 std::string h_name, h_title;
170 h_name = str(boost::format(
"ncev_Thr%1%MeV") % id1);
171 h_title = str(boost::format(
"Number of hits in event (E > %1% MeV)") % id1);
172 TH1F* h =
new TH1F(h_name.c_str(), h_title.c_str(), id2, 0, id2);
173 h->GetXaxis()->SetTitle(
"Number of hits");
177 for (
int i = 0; i < ECL_CRATES; i++) {
179 std::string h_name, h_title;
180 h_name = str(boost::format(
"time_crate_%1%_Thr1GeV") % (crate));
181 h_title = str(boost::format(
"Reconstructed time for ECL crate #%1% with E > 1 GeV") % (crate));
182 TH1F* h =
new TH1F(h_name.c_str(), h_title.c_str(), 400, -100, 100);
183 h->GetXaxis()->SetTitle(
"Time [ns]");
188 if (
id !=
"all" &&
id !=
"psd" &&
id !=
"logic" &&
id !=
"rand" &&
id !=
"dphy" &&
id !=
"other")
189 B2WARNING(
"Waveform Options are not correctly assigned. They must be 'all', 'psd', 'logic', 'rand', 'dphy', 'other'!");
191 std::string h_cell_name;
192 if (
id ==
"other") h_title =
"Unexpectedly saved waveforms";
193 if (
id ==
"psd") h_title =
"#frac{Saved}{Expected} waveforms for high-energy hits (E > 50 MeV)";
194 if (
id ==
"logic") h_title =
"#frac{Saved}{Expected} waveforms for every 1000th event";
195 if (
id ==
"rand") h_title =
"#frac{Saved}{Expected} waveforms for random trigger events";
196 if (
id ==
"dphy") h_title =
"#frac{Saved}{Expected} waveforms for delayed bhabha (DPHY) events";
197 if (
id ==
"all") h_title =
"#frac{Saved}{Expected} waveforms for all events";
198 h_cell_name = str(boost::format(
"wf_cid_%1%") % (
id));
199 TH1F* h_cell =
new TH1F(h_cell_name.c_str(), h_title.c_str(),
201 h_cell->GetXaxis()->SetTitle(
"Cell ID");
203 h_cell_psd_norm =
new TH1F(
"psd_cid",
"Normalization to psd hits for cid",
207 h_evtot_logic =
new TH1F(
"event_logic",
"Event bank for logic", 1, 0, 1);
210 h_evtot_rand =
new TH1F(
"event_rand",
"Event bank for rand", 1, 0, 1);
213 h_evtot_dphy =
new TH1F(
"event_dphy",
"Event bank for dphy", 1, 0, 1);
220 h_trigtag2_trigid =
new TH2F(
"trigtag2_trigid",
"Internal data consistency vs crate. 0-good, 1-data corruption",
221 52, 1, 53, 11, -1, 10);
228 h_pedmean_cellid->GetYaxis()->SetTitle(
"Ped. average (ADC units, #approx 0.05 MeV)");
230 h_pedrms_cellid =
new TProfile(
"pedrms_cellid",
"Pedestal stddev vs Cell ID",
233 h_pedrms_cellid->GetYaxis()->SetTitle(
"Ped. stddev (ADC units, #approx 0.05 MeV)");
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)");
240 h_trigtime_trigid =
new TH2F(
"trigtime_trigid",
"Trigger time vs Crate ID", 52, 1, 53, 145, 0, 145);
244 h_pi0_mass =
new TH1F(
"ecl_pi0_mass",
"ecl_pi0_mass", 120, 0.08, 0.20);
260 if (!
mapper.initFromDB()) B2FATAL(
"ECL DQM: Can't initialize eclChannelMapper");
284 for (TH1F* histogram :
h_cids)
286 for (TH1F* histogram :
h_edeps)
292 for (TH1F* histogram :
h_ncevs)
294 for (TH1F* histogram :
h_cells)
310 for (
auto& value :
ecltot) value = 0;
311 for (
auto& value :
nhits) value = 0;
312 bool bhatrig =
false;
316 try { bhatrig =
m_l1Trigger->testInput(
"bha_delay"); }
317 catch (
const std::exception&) { bhatrig =
false; }
334 int i = aECLDigit.getCellId() - 1;
336 if (aECLDigit.getAmp() > 2.e04 && aECLDigit.getQuality() == 3)
h_bad_quality->Fill(aECLDigit.getCellId());
337 if (aECLDigit.getAmp() >= (
v_totalthrApsd[i] / 4 * 4)) NDigits ++;
339 if (
id !=
"psd")
continue;
340 else if (
id ==
"psd" && (
m_iEvent % 1000 == 999 ||
349 double itrg = aECLTrig.getTimeTrig();
352 int tg = (int)itrg - 2 * ((
int)itrg / 8);
354 trigtag1 += aECLTrig.getTrigTag();
362 if (compar == trigtag1) flagtag = 0;
367 int cid = aECLCalDigit.getCellId();
368 double energy = aECLCalDigit.getEnergy();
369 double timing = aECLCalDigit.getTime();
373 if (energy > thrGeV) {
380 auto thrGeV = thr.value() / 1000.;
381 if (energy > thrGeV)
ecltot[thr.index()] += energy;
385 auto thrGeV = thr.value() / 1000.;
386 if (energy > thrGeV) {
387 if (cid > ECL_FWD_CHANNELS && cid <= ECL_FWD_CHANNELS + ECL_BARREL_CHANNELS)
h_time_barrels[thr.index()]->Fill(timing);
396 for (
const auto& h :
h_edeps | indexed(0)) {
397 h.value()->Fill(
ecltot[h.index()]);
400 for (
const auto& h :
h_ncevs | indexed(0)) {
401 h.value()->Fill(
nhits[h.index()]);
405 int i = aECLDsp.getCellId() - 1;
423 const auto& index = iter.index();
424 const auto& wf_opt = iter.value();
425 if (wf_opt !=
"all" && wf_opt !=
"psd" && wf_opt !=
"logic" && wf_opt !=
"rand" && wf_opt !=
"dphy" && wf_opt !=
"other")
continue;
428 else if (wf_opt ==
"logic" &&
m_iEvent % 1000 != 999)
continue;
430 else if (wf_opt ==
"dphy" && (
m_iEvent % 1000 == 999 || !bhatrig))
continue;
433 h_cells[index]->Fill(aECLDsp.getCellId());
461 const std::string trg_identifier =
"software_trigger_cut&skim&accept_hadron";
463 if (!result.isValid())
return false;
466 const std::map<std::string, int>& results = result->getResults();
467 if (results.find(trg_identifier) == results.end())
return false;
470 if (!accepted)
return false;
473 if (!particles.isValid())
return false;
475 for (
unsigned int i = 0; i < particles->getListSize(); i++) {
476 const Particle* part = particles->getParticle(i);
477 auto inv_mass = Variable::eclClusterOnlyInvariantMass(part);
double m_PedestalMean[ECLElementNumbers::c_NCrystals]
Pedestal average values.
TH1F * h_trigtag1
Histogram: Trigger tag flag #1.
std::vector< double > nhits
Container for channel multiplicity.
bool fillInvMassHistogram()
Fill entries for pi0 invariant mass distribution.
TH1F * h_bad_quality
Histogram: Cell IDs w/ bad fit quality flag.
double m_PedestalRms[ECLElementNumbers::c_NCrystals]
Pedestal rms error values.
bool m_DPHYTTYP
Flag to select events triggered by delayed bhabha.
std::vector< double > m_TimingThresholds
Parameters for timing histograms.
TH1F * h_adc_hits
Histogram: Fraction of digits above ADC threshold.
StoreArray< ECLDsp > m_ECLDsps
StoreArray ECLDsp.
virtual void initialize() override
Initialize the module.
virtual ~ECLDQMModule()
Destructor.
StoreArray< ECLCalDigit > m_ECLCalDigits
StoreArray ECLCalDigit.
StoreObjPtr< TRGSummary > m_l1Trigger
StoreObjPtr TRGSummary.
TH1F * h_cell_psd_norm
Histogram: Normalize to psd hits for CellID.
virtual void event() override
Event processor.
TH2F * h_trigtime_trigid
Histogram: Trigger time vs.
ECLDQMModule()
< derived from HistoModule class.
std::vector< TH1F * > h_edeps
Histogram vector: Total energy.
virtual void endRun() override
Call when a run ends.
std::vector< TH1F * > h_ncevs
Histogram vector: Channel multiplicity.
virtual void terminate() override
Terminate.
TH1F * h_evtot
Histogram: Total event no (auxiliary) to normalize hit map .
int m_iEvent
Global event number.
ECL::ECLGeometryPar * m_geom
Geometry.
StoreArray< ECLDigit > m_ECLDigits
StoreArray ECLDigit.
std::vector< double > ecltot
Container for energy.
std::vector< TH1F * > h_time_barrels
Histogram vector: Reconstructed time for barrel.
double m_EnergyUpperThr
Upper threshold of energy deposition in event, [GeV].
TProfile * h_pedmean_cellid
Histogram: Pedestal Average vs.
TH1F * h_evtot_rand
Histogram: Event no for rand (auxiliary) to normalize rand waveform flow.
std::vector< TH1F * > h_time_crate_Thr1GeV
Histogram vector: Reconstructed signal time for all ECL crates above the threshold = 1 GeV.
TH1F * h_time_crate_Thr1GeV_large
Histogram: Entries with crate time offsets > 100 ns (E > 1 GeV).
virtual void beginRun() override
Call when a run begins.
std::string m_histogramDirectoryName
Histogram directory in ROOT file.
ECL::ECLChannelMapper mapper
ECL channel mapper.
TH1F * h_evtot_logic
Histogram: Event no for logic (auxiliary) to normalize logic waveform flow.
std::vector< double > m_HitNumberUpperLimits
Parameters for number of hits histograms.
std::vector< double > m_TotalEnergyThresholds
Parameters for histograms w/ total energy.
TProfile * h_pedrms_thetaid
Histogram: Pedestal rms error vs.
std::vector< std::string > m_WaveformOption
Parameters for waveform histograms.
int m_DspArray[ECLElementNumbers::c_NCrystals][31]
WF sampling points for digit array.
TH1F * h_evtot_dphy
Histogram: Event no for dphy (auxiliary) to normalize dphy waveform flow.
std::vector< double > m_HitThresholds
Parameters for hit occ.
TH1F * h_quality
Histogram: Fit quality flag (0 - good, 1 - large amplitude, 3 - bad chi2).
StoreArray< ECLTrig > m_ECLTrigs
StoreArray ECLTrig.
StoreObjPtr< EventMetaData > m_eventmetadata
StoreObjPtr EventMetaData.
TH2F * h_trigtag2_trigid
Histogram: Trigger tag flag #2 vs.
std::vector< int > v_totalthrApsd
Vector to store psd wf amplitude threshold.
TH1F * h_pi0_mass
Histogram: pi0 mass.
TH1F * h_quality_other
Histogram: Fit quality flag for waveform type 'other'.
std::vector< TH1F * > h_cids
Histogram vector: Hit map.
TProfile * h_pedrms_cellid
Histogram: Pedestal rms error vs.
std::string m_pi0PListName
Name of the pi0 particle list.
DBObjPtr< ECLCrystalCalib > m_calibrationThrApsd
PSD waveform amplitude threshold.
std::vector< TH1F * > h_cells
Histogram vector: Waveforms vs CellID.
virtual void defineHisto() override
Function to define histograms.
std::vector< TH1F * > h_time_endcaps
Histogram vector: Reconstructed time for endcaps.
Class to store ECL digitized hits (output of ECLDigi) relation to ECLHit filled in ecl/modules/eclDig...
int getAmp() const
Get Fitting Amplitude.
int getQuality() const
Get Fitting Quality.
static ECLDigit * getByCellID(int cid)
Find ECLDigit by Cell ID using linear search.
static ECLGeometryPar * Instance()
Static method to get a reference to the ECLGeometryPar instance.
HistoModule()
Constructor.
void setDescription(const std::string &description)
Sets the description of the module.
void setPropertyFlags(unsigned int propertyFlags)
Sets the flags for the module properties.
@ c_ParallelProcessingCertified
This module can be run in parallel processing mode safely (All I/O must be done through the data stor...
Class to store reconstructed particles.
Type-safe access to single objects in the data store.
@ TTYP_DPHY
delayed physics events for background
@ TTYP_POIS
poisson random trigger
@ TTYP_RAND
random trigger events
static const double GeV
Standard of [energy, momentum, mass].
void addParam(const std::string &name, T ¶mVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
double sqrt(double a)
sqrt for double
@ c_accept
Accept this event.
const int c_NCrystals
Number of crystals.
Abstract base class for different kinds of events.