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
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  {
93  fCellHitLabel = pset.get< std::string >("CellHitLabel");
94  fSliceLabel = pset.get< std::string >("SliceLabel");
95  fTrackLabel = pset.get< std::string >("TrackLabel");
96  fRawDataLabel = pset.get< std::string >("RawDataLabel");
97  fTrackSyncMetric = pset.get< bool >("TrackSyncMetric");
98  fTimeLow = pset.get< double >("TimeLow");
99  fTimeHigh = pset.get< double >("TimeHigh");
100  fNoisyHitMin = pset.get< int >("NoisyHitMin");
101 
102  produces< sumdata::EventQuality >();
103  }
104 
105  //......................................................................
107  {
108  }
109 
110  //......................................................................
112  {
113  std::unique_ptr<sumdata::EventQuality> spillqual(new sumdata::EventQuality);
114 
116 
117  //get all hits info
119  evt.getByLabel(fCellHitLabel, chits);
120 
123 
125  evt.getByLabel(fRawDataLabel,daqheader);
126  if( !daqheader.failedToGet() )
127  FillRawVars(spillqual.get(),*daqheader);
128 
129  FillCountVars(spillqual.get(), *chits);
130 
131  //skip caclulating FD track sync variables
132  if(geom->DetId() == novadaq::cnv::kNEARDET){
133  evt.put(std::move(spillqual));
134  return; //skip caclulating FD track sync variables
135  }
136 
137  evt.getByLabel(fSliceLabel, slices);
138 
139  if(fTrackSyncMetric) evt.getByLabel(fTrackLabel, tracks);
140 
141  if(fTrackSyncMetric)
142  FillSyncVars(spillqual.get(), *slices, *tracks);
143  else
144  FillSyncVars(spillqual.get(), *slices, {});
145 
146 
147  //save info to event tree
148  evt.put(std::move(spillqual));
149 
150  }//end event loop
151 
152  //..............................................................
154  const std::vector<rb::CellHit>& chits) const
155  {
159 
160  unsigned int horizhits = 0; // All hits for horizontal modules
161  unsigned int ndcm3hits = 0; // DCM3 hits for horizontal modules
162 
163  unsigned int ncells_outtime = 0; // # of cells in out-beam-windomw
164 
165  // Chosen an arbitrary diblock/dcm/feb/pixel to encode a dchan in
166  // order to get diblock type (instead of hardcoding). This allows
167  // one to set the max total dcm value for either ND or FD
168  const daqchannelmap::dchan dchan = cmap->encodeDChan(geom->DetId(), 1, 1, 0, 0);
169  const daqchannelmap::DiBlock_TYPE db_t = cmap->getDiBlockType(dchan);
170 
171  // Irrespective of ND Muon Catcher, this sets max number of DCM within
172  // detector diblock to set limits on loop over DCMs later on
173  const unsigned int totdcm = cmap->getNumberOfDCM(db_t);
174  const unsigned int totdb = cmap->getTotalNumberOfDiblocks();
175 
176  std::vector<std::vector<unsigned int> > ndcmhits(totdb);
177  for(auto & itr : ndcmhits) itr.resize(totdcm, 0);
178 
179  std::map< unsigned int,
180  std::map<unsigned int, unsigned int> > dbapdmap; // Map of diblock/APD/nhits
181 
182  for(const rb::CellHit& chit: chits){
183 
184  const unsigned int dcm = cmap->getDCM(chit.DaqChannel());
185  const unsigned int db = cmap->getDiBlock(chit.DaqChannel());
186  const unsigned int apd = cmap->getFEB(chit.DaqChannel());
187  // Offset by 1 as db and dcm accounting starts at 1
188  ++ndcmhits[db-1][dcm-1];
189 
190  // Only look at hits outside of beam window
191  if(chit.TNS() < fTimeLow || chit.TNS() > fTimeHigh ){
192  ++ncells_outtime; // total number of cells out of time
193 
194  // Cells 95, 87, 79 and 71 in Y-view in ND are impacted by
195  // lights on the catwalk. See K. Sachdev doc-db 12426
196  if(geom->DetId() == novadaq::cnv::kNEARDET &&
197  chit.Plane() < geom->FirstPlaneInMuonCatcher() &&
198  chit.View() == geo::kY){
199  ++horizhits;
200  int cell = chit.Cell();
201  if(cell == 95 || cell == 87 || cell == 79 || cell == 71)
202  ++ndcm3hits;
203  } // DB horizontal plane only, excluding Muon Catcher
204 
205  // count number of hits for each apd into ndbapd map
206  ++dbapdmap[db][apd];
207  } // outside beam window
208  } // end loop over all hits
209 
210  unsigned int ndeaddcms = 0;
211 
212  // Loop over di-blocks
213  for(unsigned int dbIdx = 0; dbIdx < totdb; ++dbIdx){
214  if(geom->DetId() == novadaq::cnv::kNEARDET && (dbIdx < (totdb-1))){
215  for(unsigned int dcmIdx = 0; dcmIdx < totdcm; ++dcmIdx){
216  if(ndcmhits[dbIdx][dcmIdx] == 0) ndeaddcms++;
217  }
218  } // end loop for ND minus Muon Catcher
219  else if(geom->DetId() != novadaq::cnv::kNEARDET){
220  for(unsigned int dcmIdx = 0; dcmIdx < totdcm; ++dcmIdx){
221  if(ndcmhits[dbIdx][dcmIdx] == 0) ++ndeaddcms;
222  }
223  } // end loop for NOT ND
224  } // end loop over diblocks
225  // Special case for the 2 Near Det Muon Catcher DCMS
226  if(geom->DetId() == novadaq::cnv::kNEARDET){
227  for(unsigned int dcmIdx = 0; dcmIdx < 2; ++dcmIdx){
228  if(ndcmhits[totdb-1][dcmIdx] == 0) ++ndeaddcms;
229  } // end loop for ND Muon Catcher DCMs
230  } // end loop for ND
231 
232  // Calculate fraction of partial DCM3 hits in horizontal modules to
233  // detect lights-on effect in ND
234  float filter_lightson = 0.; // fraction of partial DCM3 hits in horizontal modules
235  if(horizhits > 0) filter_lightson = float(ndcm3hits)/horizhits;
236 
237  // Noisy APD if it has at least fNoisyHitMin hits.
238  // Original value of fNoisyHitMin = 28 for ND
239  // See X. Bu doc-db 12316
240  // Total number of noisy apd in ND
241  int napd_hitmin_tot = 0;
242 
243  if(geom->DetId() == novadaq::cnv::kNEARDET){
244  for(const auto dbItr : dbapdmap){
245  for(const auto apdItr : dbItr.second){
246  if(apdItr.second > fNoisyHitMin) napd_hitmin_tot++;
247  } // end loop APD/Value pair
248  } // end loop dbapdmap
249  } //end detector check
250 
251  //save the values for dq flags
252  spillqual->nmissingdcmslg = livegeom->NDropouts();
253  spillqual->nmissingdcms = ndeaddcms;
254  spillqual->fracdcm3hits = filter_lightson;
255  spillqual->nouttimehits = ncells_outtime;
256  spillqual->nnoisyapds = napd_hitmin_tot;
257  }
258 
259  //..............................................................
261  const std::vector<rb::Cluster>& slices,
262  const std::vector<rb::Track>& tracks) const
263  {
264  // Which DCMs are populated at all, so we don't penalize at edges and for
265  // missing DCMs (that's a different cut)
266  std::set<DCMId> pop;
267 
268  for(const rb::Cluster& slice: slices) FillPopulatedDCMs(slice, pop);
269 
270  // How many edge crossings we considered
271  int nedge = 0;
272  // How many had a match
273  int nmatch = 0;
274 
275  if(fTrackSyncMetric){
276  for(const rb::Track& track: tracks)
277  AccumulateSyncMetric(track, pop, nedge, nmatch);
278  }
279  else{
280  for(const rb::Cluster& slice: slices){
281  if(!slice.IsNoise()) AccumulateSyncMetric(slice, pop, nedge, nmatch);
282  }
283  }
284 
285  spillqual->dcmedgematchfrac = (nedge ? nmatch/double(nedge) : -1);
286  }
287 
289 
290  //..............................................................
292  {
293  switch(d){
294  case kLeft: return kRight;
295  case kRight: return kLeft;
296  case kUp: return kDown;
297  case kDown: return kUp;
298  default: assert(0 && "Not reached");
299  }
300  }
301 
302  //..............................................................
304  {
305  switch(d){
306  case kLeft: return DCMId(dcm.view, dcm.x-1, dcm.y );
307  case kRight: return DCMId(dcm.view, dcm.x+1, dcm.y );
308  case kUp: return DCMId(dcm.view, dcm.x, dcm.y+1);
309  case kDown: return DCMId(dcm.view, dcm.x, dcm.y-1);
310  default: assert(0 && "Not reached");
311  }
312  }
313 
314  struct Edges
315  {
316  Edges() : flag() {}
317  bool flag[4];
318  };
319 
320  //..............................................................
322  std::set<DCMId>& pop) const
323  {
324  for(unsigned int chitIdx = 0; chitIdx < clust.NCell(); ++chitIdx){
325  art::Ptr<rb::CellHit> chit = clust.Cell(chitIdx);
326 
327  const int dcmx = chit->Plane() / 64;
328  const int dcmy = chit->Cell() / 64;
329 
330  pop.insert(DCMId(chit->View(), dcmx, dcmy));
331  }
332  }
333 
334  //..............................................................
336  const std::set<DCMId>& pop,
337  int& nedge, int& nmatch) const
338  {
339  std::map<DCMId, Edges> edges;
340 
341  for(unsigned int chitIdx = 0; chitIdx < clust.NCell(); ++chitIdx){
342  art::Ptr<rb::CellHit> chit = clust.Cell(chitIdx);
343 
344  const int modplane = chit->Plane() % 64;
345  const int modcell = chit->Cell() % 64;
346 
347  const int dcmx = chit->Plane() / 64;
348  const int dcmy = chit->Cell() / 64;
349 
350  const DCMId key(chit->View(), dcmx, dcmy);
351 
352  if(modplane == 0 || modplane == 1) edges[key].flag[kLeft] = true;
353  if(modplane == 62 || modplane == 63) edges[key].flag[kRight] = true;
354  if(modcell == 0) edges[key].flag[kDown] = true;
355  if(modcell == 63) edges[key].flag[kUp] = true;
356  } // end for chitIdx
357 
358  // Make a copy because [] access to edges in the loop below might mutate
359  // what we're looping over?
360  std::map<DCMId, Edges> edges_it = edges;
361 
362  for(auto it: edges_it){
363  const DCMId key = it.first;
364  const Edges val = it.second;
365 
366  for(Dir dir: {kLeft, kRight, kUp, kDown}){
367  if(!val.flag[dir]) continue;
368 
369  const DCMId nei = Move(key, dir);
370  // Only count if the neighbouring DCM has hits at all
371  if(!pop.count(nei)) continue;
372 
373  if(edges[nei].flag[Opposite(dir)]){
374  ++nedge;
375  ++nmatch;
376  }
377  else{
378  // If there's no match, count two chances, because a good match
379  // will have been counted from both sides.
380  nedge += 2;
381  }
382  } // end for dir
383  } // end for it
384  }
385 
386  //-----------------------------------------------------
388  const rawdata::DAQHeader& daqheader)
389  {
390  spillqual->nmicroslices = daqheader.TotalMicroSlices();
391  }//end for FillRawVars
392 
393 } // end namespace dqsf
394 ////////////////////////////////////////////////////////////////////////
395 
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
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
unsigned short Cell() const
Definition: CellHit.h:40
void FillPopulatedDCMs(const rb::Cluster &clust, std::set< DCMId > &pop) const
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
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
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
TDirectory * dir
Definition: macro.C:5
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
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.
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:196
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.