DQSpillFlags_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 ///
3 /// \brief Calculate the spill level data quality flags and
4 /// save to sumdata::EventQuality
5 /// \authors xbbu@fnal.gov
6 /// \date 11/10/2014
7 ///
8 ////////////////////////////////////////////////////////////////////////
9 
10 //C++ includes
11 #include <cmath>
12 #include <map>
13 #include <string>
14 #include <vector>
15 
16 // Framework includes
21 #include "art_root_io/TFileDirectory.h"
22 #include "art_root_io/TFileService.h"
24 #include "fhiclcpp/ParameterSet.h"
25 
26 // NOvASoft includes
27 #include "Calibrator/Calibrator.h"
29 #include "DAQChannelMap/DAQChannelMap.h"
30 #include "Geometry/Geometry.h"
31 #include "Geometry/LiveGeometry.h"
32 #include "NovaDAQConventions/DAQConventions.h"
33 #include "RecoBase/CellHit.h"
34 #include "RecoBase/Cluster.h"
35 #include "RecoBase/Track.h"
37 #include "RawData/DAQHeader.h"
38 
39 namespace dqsf
40 {
41  struct DCMId
42  {
43  DCMId(geo::View_t v, int _x, int _y) : view(v), x(_x), y(_y) {}
45  int x, y;
46 
47  bool operator<(const DCMId rhs) const
48  {
49  return std::make_tuple(view, x, y) < std::make_tuple(rhs.view, rhs.x, rhs.y);
50  }
51  };
52 
53  class DQSpillFlags : public art::EDProducer
54  {
55  public:
56  explicit DQSpillFlags(fhicl::ParameterSet const& pset);
57  virtual ~DQSpillFlags();
58 
59  void produce(art::Event& evt);
60 
61  protected:
62  void FillCountVars(sumdata::EventQuality* spillqual,
63  const std::vector<rb::CellHit>& chits) const;
64  void FillSyncVars(sumdata::EventQuality* spillqual,
65  const std::vector<rb::Cluster>& slices,
66  const std::vector<rb::Track>& tracks) const;
67  void FillRawVars(sumdata::EventQuality* spillqual,
68  const rawdata::DAQHeader& daqheader);
69 
70  void FillPopulatedDCMs(const rb::Cluster& clust,
71  std::set<DCMId>& pop) const;
72 
73  void AccumulateSyncMetric(const rb::Cluster& clust,
74  const std::set<DCMId>& pop,
75  int& nedge, int& nmatch) const;
76 
77  std::string fCellHitLabel; ///< Label for cell hits
78  std::string fSliceLabel; ///< Label for slices
79  std::string fTrackLabel; ///< Label for tracks
80  std::string fRawDataLabel; ///< Label for raw data
81 
82  bool fTrackSyncMetric; ///< Use track hits rather than all hits
83 
84  double fTimeLow; ///< Hit time of beam trigger window low end
85  double fTimeHigh; ///< Hit time of beam trigger window high end
86  unsigned int fNoisyHitMin; ///< Miniumum number of hits to be considered a noisy APD
87  };
88 
89 
90  //..............................................................
92  : EDProducer(pset)
93  {
94  fCellHitLabel = pset.get< std::string >("CellHitLabel");
95  fSliceLabel = pset.get< std::string >("SliceLabel");
96  fTrackLabel = pset.get< std::string >("TrackLabel");
97  fRawDataLabel = pset.get< std::string >("RawDataLabel");
98  fTrackSyncMetric = pset.get< bool >("TrackSyncMetric");
99  fTimeLow = pset.get< double >("TimeLow");
100  fTimeHigh = pset.get< double >("TimeHigh");
101  fNoisyHitMin = pset.get< int >("NoisyHitMin");
102 
103  produces< sumdata::EventQuality >();
104  }
105 
106  //......................................................................
108  {
109  }
110 
111  //......................................................................
113  {
114  std::unique_ptr<sumdata::EventQuality> spillqual(new sumdata::EventQuality);
115 
117 
118  //get all hits info
120  evt.getByLabel(fCellHitLabel, chits);
121 
124 
126  evt.getByLabel(fRawDataLabel,daqheader);
127  if( !daqheader.failedToGet() )
128  FillRawVars(spillqual.get(),*daqheader);
129 
130  FillCountVars(spillqual.get(), *chits);
131 
132  //skip caclulating FD track sync variables
133  if(geom->DetId() == novadaq::cnv::kNEARDET){
134  evt.put(std::move(spillqual));
135  return; //skip caclulating FD track sync variables
136  }
137 
138  evt.getByLabel(fSliceLabel, slices);
139 
140  if(fTrackSyncMetric) evt.getByLabel(fTrackLabel, tracks);
141 
142  if(fTrackSyncMetric)
143  FillSyncVars(spillqual.get(), *slices, *tracks);
144  else
145  FillSyncVars(spillqual.get(), *slices, {});
146 
147 
148  //save info to event tree
149  evt.put(std::move(spillqual));
150 
151  }//end event loop
152 
153  //..............................................................
155  const std::vector<rb::CellHit>& chits) const
156  {
160 
161  unsigned int horizhits = 0; // All hits for horizontal modules
162  unsigned int ndcm3hits = 0; // DCM3 hits for horizontal modules
163 
164  unsigned int ncells_outtime = 0; // # of cells in out-beam-windomw
165 
166  // Chosen an arbitrary diblock/dcm/feb/pixel to encode a dchan in
167  // order to get diblock type (instead of hardcoding). This allows
168  // one to set the max total dcm value for either ND or FD
169  const daqchannelmap::dchan dchan = cmap->encodeDChan(geom->DetId(), 1, 1, 0, 0);
170  const daqchannelmap::DiBlock_TYPE db_t = cmap->getDiBlockType(dchan);
171 
172  // Irrespective of ND Muon Catcher, this sets max number of DCM within
173  // detector diblock to set limits on loop over DCMs later on
174  const unsigned int totdcm = cmap->getNumberOfDCM(db_t);
175  const unsigned int totdb = cmap->getTotalNumberOfDiblocks();
176 
177  std::vector<std::vector<unsigned int> > ndcmhits(totdb);
178  for(auto & itr : ndcmhits) itr.resize(totdcm, 0);
179 
180  std::map< unsigned int,
181  std::map<unsigned int, unsigned int> > dbapdmap; // Map of diblock/APD/nhits
182 
183  for(const rb::CellHit& chit: chits){
184 
185  const unsigned int dcm = cmap->getDCM(chit.DaqChannel());
186  const unsigned int db = cmap->getDiBlock(chit.DaqChannel());
187  const unsigned int apd = cmap->getFEB(chit.DaqChannel());
188  // Offset by 1 as db and dcm accounting starts at 1
189  ++ndcmhits[db-1][dcm-1];
190 
191  // Only look at hits outside of beam window
192  if(chit.TNS() < fTimeLow || chit.TNS() > fTimeHigh ){
193  ++ncells_outtime; // total number of cells out of time
194 
195  // Cells 95, 87, 79 and 71 in Y-view in ND are impacted by
196  // lights on the catwalk. See K. Sachdev doc-db 12426
197  if(geom->DetId() == novadaq::cnv::kNEARDET &&
198  chit.Plane() < geom->FirstPlaneInMuonCatcher() &&
199  chit.View() == geo::kY){
200  ++horizhits;
201  int cell = chit.Cell();
202  if(cell == 95 || cell == 87 || cell == 79 || cell == 71)
203  ++ndcm3hits;
204  } // DB horizontal plane only, excluding Muon Catcher
205 
206  // count number of hits for each apd into ndbapd map
207  ++dbapdmap[db][apd];
208  } // outside beam window
209  } // end loop over all hits
210 
211  unsigned int ndeaddcms = 0;
212 
213  // Loop over di-blocks
214  for(unsigned int dbIdx = 0; dbIdx < totdb; ++dbIdx){
215  if(geom->DetId() == novadaq::cnv::kNEARDET && (dbIdx < (totdb-1))){
216  for(unsigned int dcmIdx = 0; dcmIdx < totdcm; ++dcmIdx){
217  if(ndcmhits[dbIdx][dcmIdx] == 0) ndeaddcms++;
218  }
219  } // end loop for ND minus Muon Catcher
220  else if(geom->DetId() != novadaq::cnv::kNEARDET){
221  for(unsigned int dcmIdx = 0; dcmIdx < totdcm; ++dcmIdx){
222  if(ndcmhits[dbIdx][dcmIdx] == 0) ++ndeaddcms;
223  }
224  } // end loop for NOT ND
225  } // end loop over diblocks
226  // Special case for the 2 Near Det Muon Catcher DCMS
227  if(geom->DetId() == novadaq::cnv::kNEARDET){
228  for(unsigned int dcmIdx = 0; dcmIdx < 2; ++dcmIdx){
229  if(ndcmhits[totdb-1][dcmIdx] == 0) ++ndeaddcms;
230  } // end loop for ND Muon Catcher DCMs
231  } // end loop for ND
232 
233  // Calculate fraction of partial DCM3 hits in horizontal modules to
234  // detect lights-on effect in ND
235  float filter_lightson = 0.; // fraction of partial DCM3 hits in horizontal modules
236  if(horizhits > 0) filter_lightson = float(ndcm3hits)/horizhits;
237 
238  // Noisy APD if it has at least fNoisyHitMin hits.
239  // Original value of fNoisyHitMin = 28 for ND
240  // See X. Bu doc-db 12316
241  // Total number of noisy apd in ND
242  int napd_hitmin_tot = 0;
243 
244  if(geom->DetId() == novadaq::cnv::kNEARDET){
245  for(const auto dbItr : dbapdmap){
246  for(const auto apdItr : dbItr.second){
247  if(apdItr.second > fNoisyHitMin) napd_hitmin_tot++;
248  } // end loop APD/Value pair
249  } // end loop dbapdmap
250  } //end detector check
251 
252  //save the values for dq flags
253  spillqual->nmissingdcmslg = livegeom->NDropouts();
254  spillqual->nmissingdcms = ndeaddcms;
255  spillqual->fracdcm3hits = filter_lightson;
256  spillqual->nouttimehits = ncells_outtime;
257  spillqual->nnoisyapds = napd_hitmin_tot;
258  }
259 
260  //..............................................................
262  const std::vector<rb::Cluster>& slices,
263  const std::vector<rb::Track>& tracks) const
264  {
265  // Which DCMs are populated at all, so we don't penalize at edges and for
266  // missing DCMs (that's a different cut)
267  std::set<DCMId> pop;
268 
269  for(const rb::Cluster& slice: slices) FillPopulatedDCMs(slice, pop);
270 
271  // How many edge crossings we considered
272  int nedge = 0;
273  // How many had a match
274  int nmatch = 0;
275 
276  if(fTrackSyncMetric){
277  for(const rb::Track& track: tracks)
278  AccumulateSyncMetric(track, pop, nedge, nmatch);
279  }
280  else{
281  for(const rb::Cluster& slice: slices){
282  if(!slice.IsNoise()) AccumulateSyncMetric(slice, pop, nedge, nmatch);
283  }
284  }
285 
286  spillqual->dcmedgematchfrac = (nedge ? nmatch/double(nedge) : -1);
287  }
288 
290 
291  //..............................................................
293  {
294  switch(d){
295  case kLeft: return kRight;
296  case kRight: return kLeft;
297  case kUp: return kDown;
298  case kDown: return kUp;
299  default: assert(0 && "Not reached");
300  }
301  }
302 
303  //..............................................................
305  {
306  switch(d){
307  case kLeft: return DCMId(dcm.view, dcm.x-1, dcm.y );
308  case kRight: return DCMId(dcm.view, dcm.x+1, dcm.y );
309  case kUp: return DCMId(dcm.view, dcm.x, dcm.y+1);
310  case kDown: return DCMId(dcm.view, dcm.x, dcm.y-1);
311  default: assert(0 && "Not reached");
312  }
313  }
314 
315  struct Edges
316  {
317  Edges() : flag() {}
318  bool flag[4];
319  };
320 
321  //..............................................................
323  std::set<DCMId>& pop) const
324  {
325  for(unsigned int chitIdx = 0; chitIdx < clust.NCell(); ++chitIdx){
326  art::Ptr<rb::CellHit> chit = clust.Cell(chitIdx);
327 
328  const int dcmx = chit->Plane() / 64;
329  const int dcmy = chit->Cell() / 64;
330 
331  pop.insert(DCMId(chit->View(), dcmx, dcmy));
332  }
333  }
334 
335  //..............................................................
337  const std::set<DCMId>& pop,
338  int& nedge, int& nmatch) const
339  {
340  std::map<DCMId, Edges> edges;
341 
342  for(unsigned int chitIdx = 0; chitIdx < clust.NCell(); ++chitIdx){
343  art::Ptr<rb::CellHit> chit = clust.Cell(chitIdx);
344 
345  const int modplane = chit->Plane() % 64;
346  const int modcell = chit->Cell() % 64;
347 
348  const int dcmx = chit->Plane() / 64;
349  const int dcmy = chit->Cell() / 64;
350 
351  const DCMId key(chit->View(), dcmx, dcmy);
352 
353  if(modplane == 0 || modplane == 1) edges[key].flag[kLeft] = true;
354  if(modplane == 62 || modplane == 63) edges[key].flag[kRight] = true;
355  if(modcell == 0) edges[key].flag[kDown] = true;
356  if(modcell == 63) edges[key].flag[kUp] = true;
357  } // end for chitIdx
358 
359  // Make a copy because [] access to edges in the loop below might mutate
360  // what we're looping over?
361  std::map<DCMId, Edges> edges_it = edges;
362 
363  for(auto it: edges_it){
364  const DCMId key = it.first;
365  const Edges val = it.second;
366 
367  for(Dir dir: {kLeft, kRight, kUp, kDown}){
368  if(!val.flag[dir]) continue;
369 
370  const DCMId nei = Move(key, dir);
371  // Only count if the neighbouring DCM has hits at all
372  if(!pop.count(nei)) continue;
373 
374  if(edges[nei].flag[Opposite(dir)]){
375  ++nedge;
376  ++nmatch;
377  }
378  else{
379  // If there's no match, count two chances, because a good match
380  // will have been counted from both sides.
381  nedge += 2;
382  }
383  } // end for dir
384  } // end for it
385  }
386 
387  //-----------------------------------------------------
389  const rawdata::DAQHeader& daqheader)
390  {
391  spillqual->nmicroslices = daqheader.TotalMicroSlices();
392  }//end for FillRawVars
393 
394 } // end namespace dqsf
395 ////////////////////////////////////////////////////////////////////////
396 
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
int nouttimehits
of out-beam-window hits
Definition: EventQuality.h:25
virtual unsigned int getNumberOfDCM(DiBlock_TYPE dbt) const =0
How many DCMs does this diblock have?
set< int >::iterator it
std::string fTrackLabel
Label for tracks.
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
unsigned short Plane() const
Definition: CellHit.h:39
std::string fSliceLabel
Label for slices.
geo::View_t View() const
Definition: CellHit.h:41
std::string fRawDataLabel
Label for raw data.
Definition: event.h:19
A collection of associated CellHits.
Definition: Cluster.h:47
double fTimeHigh
Hit time of beam trigger window high end.
DCMId(geo::View_t v, int _x, int _y)
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
DQSpillFlags(fhicl::ParameterSet const &pset)
DEFINE_ART_MODULE(TestTMapFile)
int TotalMicroSlices() const
Definition: DAQHeader.h:27
int nmissingdcms
of missing DCMs
Definition: EventQuality.h:23
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
int nmicroslices
of micro slices
Definition: EventQuality.h:27
unsigned short Cell() const
Definition: CellHit.h:40
void FillPopulatedDCMs(const rb::Cluster &clust, std::set< DCMId > &pop) const
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:446
void FillRawVars(sumdata::EventQuality *spillqual, const rawdata::DAQHeader &daqheader)
double fracdcm3hits
fraction of DCM3 hits in horizontal modules
Definition: EventQuality.h:24
float dcmedgematchfrac
Low values mean out-of-sync detector.
Definition: EventQuality.h:29
T get(std::string const &key) const
Definition: ParameterSet.h:231
virtual unsigned int getTotalNumberOfDiblocks() const =0
How many diblocks does the detector have?
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
int evt
geo::View_t view
Near Detector in the NuMI cavern.
Float_t d
Definition: plot.C:236
bool operator<(const DCMId rhs) const
void FillCountVars(sumdata::EventQuality *spillqual, const std::vector< rb::CellHit > &chits) const
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
TDirectory * dir
Definition: macro.C:5
void geom(int which=0)
Definition: geom.C:163
void AccumulateSyncMetric(const rb::Cluster &clust, const std::set< DCMId > &pop, int &nedge, int &nmatch) const
DCMId Move(DCMId dcm, Dir d)
assert(nhit_max >=nhit_nbins)
cmap::CMap class source code
Definition: CMap.cxx:17
dcm_id_t getDCM(dchan daqchan) const
Decode the dcm ID from a dchan.
dchan encodeDChan(int detID, diblock_t diblock, dcm_id_t dcm, feb_t feb, pixel_t pixel) const
DiBlock_TYPE
Types of Diblock.
Calculate the spill level data quality flags and save to sumdata::EventQuality.
uint32_t dchan
< DAQ Channel Map Package
bool fTrackSyncMetric
Use track hits rather than all hits.
ProductID put(std::unique_ptr< PROD > &&edp, FullSemantic< Level::Run > const semantic)
Definition: DataViewImpl.h:730
unsigned int fNoisyHitMin
Miniumum number of hits to be considered a noisy APD.
double fTimeLow
Hit time of beam trigger window low end.
int nnoisyapds
of noisy APDs
Definition: EventQuality.h:26
diblock_t getDiBlock(dchan daqchan) const
Decode the diblock ID from a dchan.
virtual DiBlock_TYPE getDiBlockType(dchan chan) const =0
What format is the diblock? Only relevant in NDOS.
Encapsulate the geometry of one entire detector (near, far, ndos)
std::string fCellHitLabel
Label for cell hits.
bool failedToGet() const
Definition: Handle.h:190
void FillSyncVars(sumdata::EventQuality *spillqual, const std::vector< rb::Cluster > &slices, const std::vector< rb::Track > &tracks) const
void produce(art::Event &evt)
Dir Opposite(Dir d)
feb_t getFEB(dchan daqchan) const
Decode the feb id from a dchan.
const unsigned int FirstPlaneInMuonCatcher() const
Returns the index of the first plane contained in the muon catcher.
enum BeamMode string