DataPrepByApaModule_module.cc
Go to the documentation of this file.
1 // DataPrepByApaModule_module.cc
2 
3 // David Adams
4 // October 2019
5 //
6 // Module that obtain RawDigit and TimeStamp containers from a decoder tool, builds
7 // AdcChannelData containers and runs RawDigitPrepService on those.
8 // It may write all the decoded digits and/or time stamps and selected wires.
9 //
10 // Configuration parameters:
11 // LogLevel - Usual logging level.
12 // DecoderTool - Name of the tool used to fetch the raw digits.
13 // OutputDigitName - Name for the output digit container. If blank, digits are not written.
14 // OutputTimeStampName - Name for the output wire container. If blank, time stamps are not written.
15 // OutputWireName - Name for the output wire container. If blank, wires are not written.
16 // ChannelGroups - Process channels appearing in these groups or ranges.
17 // The group for each name is obtained from the tool channelGroups.
18 // If that fails, a range with that name is obtained from the tool
19 // channelRanges.
20 // If the entry "all" appears, then all channels are decoded at once and all
21 // channels are processed together.
22 // Otherwise each APA is read in and processed individually.
23 // The group "apas" should be defined to include all APAs.
24 // SkipEmptyChannels - If true, empty chanels (raw data length zero) are ignored
25 // DeltaTickCount - 0 - Keep channels with tick count equal to Nprim (primary tick count)
26 // > 0 - Keep channels within delta*nPrim of Nprim
27 // < 0 - No cut on tick count.
28 // BeamEventLabel - Label for the BeamEvent (trigger) data product. If blank, it is not read.
29 // ApaChannelCounts - Number of channels in each APA. Last value is used for subsequent APAs.
30 // OnlineChannelMapTool - Name of the tool that converts offline to online channel number.
31 
49 #include "TTimeStamp.h"
52 #include <iomanip>
53 #include <set>
54 #include <sstream>
55 
56 using std::cout;
57 using std::endl;
58 using std::string;
59 using std::vector;
60 using std::move;
61 using std::setw;
62 using art::ServiceHandle;
63 using art::Timestamp;
64 using raw::RDStatus;
65 using recob::Wire;
66 using std::istringstream;
67 using std::ostringstream;
68 
69 //**********************************************************************
70 
72 
73 public:
74 
75  using Index = unsigned int;
76  using Name = std::string;
77  using NameVector = std::vector<Name>;
78  using ApaNamesMap = std::map<int, NameVector>;
79  using ApaChannelSet = std::set<Index>;
80  using ApaChannelSetMap = std::map<int, ApaChannelSet>;
81 
82  // Ctor.
83  explicit DataPrepByApaModule(fhicl::ParameterSet const& pset);
84 
85  // Dtor.
87 
88  // Producer methods.
89  void produce(art::Event& evt);
90  void beginJob();
91  void endJob();
92  void reconfigure(fhicl::ParameterSet const& p);
93 
94 private:
95 
96  // Configuration parameters.
98  Name m_DecoderTool; // Name for the decoder tool
99  Name m_OutputTimeStampName; // Label for the output RDTimeStamp conainer
100  Name m_OutputDigitName; // Label for the output raw::RawDigit conainer
101  Name m_OutputWireName; ///< Second field in full label for the output wire container.
109 
110  // Split label into producer and name: PRODUCER or PRODUCER:NAME
113 
114  // Accessed services.
117 
118  // Tools.
120  std::unique_ptr<PDSPTPCDataInterfaceParent> m_pDecoderTool;
121  std::unique_ptr<IndexMapTool> m_onlineChannelMapTool;
122 
123  // Channel range names and channel sets indexed by APA.
126 
127  // Maximum size for the output containers.
131 
132  // Processed event count.
134 
135  // Skipped event count.
137 
138 };
139 
141 
142 //**********************************************************************
143 
145  const Name myname = "DataPrepByApaModule::ctor: ";
146  this->reconfigure(pset);
147  if ( m_OutputTimeStampName.size() ) {
148  if ( m_LogLevel > 0 ) {
149  cout << myname << "Module will produce RDTimeStamps with name " << m_OutputTimeStampName << endl;
150  }
151  produces<std::vector<raw::RDTimeStamp>>(m_OutputTimeStampName);
152  } else if ( m_LogLevel > 0 ) {
153  cout << myname << "Module will not produce RDTimeStamps." << endl;
154  }
155  if ( m_OutputDigitName.size() ) {
156  if ( m_LogLevel > 0 ) {
157  cout << myname << "Module will produce digits with name " << m_OutputDigitName << endl;
158  }
159  produces<std::vector<raw::RawDigit>>(m_OutputDigitName);
160  } else if ( m_LogLevel > 0 ) {
161  cout << myname << "Module will not produce RawDigits." << endl;
162  }
163  if ( m_OutputWireName.size() ) {
164  if ( m_LogLevel > 0 ) {
165  cout << myname << "Module will produce Wires with name " << m_OutputWireName << endl;
166  }
167  produces<std::vector<recob::Wire>>(m_OutputWireName);
168  } else if ( m_LogLevel > 0 ) {
169  cout << myname << "Module will not produce Wires." << endl;
170  }
171 
172  if ( m_OutputDigitName.size() ) {
173  if ( m_LogLevel > 0 ) {
174  cout << myname << "Module will produce statusess with name " << m_OutputDigitName << endl;
175  }
176  produces<std::vector<raw::RDStatus>>(m_OutputDigitName);
177  }
178 }
179 
180 //**********************************************************************
181 
183 
184 //**********************************************************************
185 
187  const string myname = "DataPrepByApaModule::reconfigure: ";
188  m_LogLevel = pset.get<int>("LogLevel");
189  m_DecoderTool = pset.get<Name>("DecoderTool");
190  m_OutputTimeStampName = pset.get<Name>("OutputTimeStampName");
191  m_OutputDigitName = pset.get<Name>("OutputDigitName");
192  m_OutputWireName = pset.get<Name>("OutputWireName", "");
193  m_ChannelGroups = pset.get<NameVector>("ChannelGroups");
194  m_BeamEventLabel = pset.get<string>("BeamEventLabel");
195  m_KeepChannels = pset.get<AdcChannelVector>("KeepChannels");
196  m_SkipChannels = pset.get<AdcChannelVector>("SkipChannels");
197  m_SkipEmptyChannels = pset.get<bool>("SkipEmptyChannels");
198  m_DeltaTickCount = pset.get<float>("DeltaTickCount");
199  m_ApaChannelCounts = pset.get<AdcChannelVector>("ApaChannelCounts");
200  pset.get_if_present<std::string>("OnlineChannelMapTool", m_OnlineChannelMapTool);
201 
202  // Display configuration.
203  if ( m_LogLevel >= 1 ) {
204  cout << myname << " LogLevel: " << m_LogLevel << endl;
205  cout << myname << " DecoderTool: " << m_DecoderTool << endl;
206  cout << myname << " OutputTimeStampName: " << m_OutputTimeStampName << endl;
207  cout << myname << " OutputDigitName: " << m_OutputDigitName << endl;
208  cout << myname << " OutputWireName: " << m_OutputWireName << endl;
209  cout << myname << " ChannelGroups: [";
210  bool first = true;
211  for ( Name rnam : m_ChannelGroups ) {
212  if ( first ) first = false;
213  else cout << ", ";
214  cout << rnam;
215  }
216  cout << "]" << endl;
217  cout << myname << " BeamEventLabel: " << m_BeamEventLabel << endl;
218  cout << myname << " OnlineChannelMapTool: " << m_OnlineChannelMapTool << endl;
219  cout << myname << " KeepChannels: " << "[";
220  first = true;
221  for ( AdcChannel ich : m_KeepChannels ) {
222  if ( first ) first = false;
223  else cout << ", ";
224  cout << ich;
225  }
226  cout << "]" << endl;
227  cout << myname << " SkipChannels: " << "[";
228  first = true;
229  for ( AdcChannel ich : m_SkipChannels ) {
230  if ( first ) first = false;
231  else cout << ", ";
232  cout << ich;
233  }
234  cout << "]" << endl;
235  cout << myname << " SkipEmptyChannels: " << (m_SkipEmptyChannels ? "true" : "false") << endl;
236  cout << myname << " DeltaTickCount: " << m_DeltaTickCount << endl;
237  cout << myname << " ApaChannelCounts: " << "[";
238  first = true;
239  for ( AdcChannel ich : m_ApaChannelCounts ) {
240  if ( first ) first = false;
241  else cout << ", ";
242  cout << ich;
243  }
244  cout << "]" << endl;
245  }
246 
247  // Fetch services.
249  if ( m_pChannelStatusProvider == nullptr ) {
250  cout << myname << "WARNING: Channel status provider not found." << endl;
251  }
253 
254  // Fetch tools.
256  if ( ptm == nullptr ) {
257  cout << myname << "ERROR: Unable to retrieve tool manager." << endl;
258  }
259  if ( m_DecoderTool.size() ) {
261  if ( ! m_pDecoderTool ) {
262  cout << myname << "ERROR: Decoder tool not found: " << m_DecoderTool << endl;
263  return;
264  }
265  }
266  if ( m_OnlineChannelMapTool.size() ) {
268  }
269  const IndexRangeTool* pcrt = ptm->getShared<IndexRangeTool>("channelRanges");
270  if ( pcrt == nullptr ) {
271  cout << myname << "ERROR: Channel range tool channelRanges not found." << endl;
272  return;
273  }
274  const IndexRangeGroupTool* pcgt = ptm->getShared<IndexRangeGroupTool>("channelGroups");
275  if ( pcgt == nullptr ) {
276  cout << myname << "WARNING: Channel group tool channelGroups not found." << endl;
277  }
278 
279  // Build the vector of channel range names from the vector of channel group names.
280  // If a group name is valid, add that group's ranges.
281  // Otherwise, assume the group name is a range name.
282  bool useGroups = pcgt != nullptr;
283  NameVector rangeNames;
284  for ( Name groupName : m_ChannelGroups ) {
285  IndexRangeGroup grp;
286  if ( useGroups ) grp = pcgt->get(groupName);
287  if ( grp.isValid() ) {
288  for ( const IndexRange& ran : grp.ranges ) {
289  rangeNames.push_back(ran.name);
290  }
291  } else {
292  rangeNames.push_back(groupName);
293  }
294  }
295 
296  // Copy the channels to skip to a set for fast lookup.
297  ApaChannelSet skipChans;
298  for ( Index icha : m_SkipChannels ) skipChans.insert(icha);
299  ApaChannelSet keepChans;
300  for ( Index icha : m_KeepChannels ) keepChans.insert(icha);
301 
302  // Sort channel ranges by APA and record the channels for each..
303  int ncrn = 0;
304  AdcChannel nchaAll = 0;
305  for ( Name crn : rangeNames ) {
306  // If any range is named "all", we keep it only.
307  int iapa = -99;
308  if ( crn == "all" ) {
309  m_apacrns.clear();
310  m_apachsets.clear();
311  iapa = -1;
312  ncrn = -1;
313  // Temporary hack for VD cold box.
314  } else if ( crn.substr(0,2) == "cr" ) {
315  iapa = 1;
316  } else {
317  string::size_type ipos = string::npos;
318  if ( crn.substr(0,3) == "apa" ) ipos = 3;
319  if ( crn.substr(0,4) == "femb" ) ipos = 4;
320  if ( ipos == string::npos ) {
321  cout << myname << "WARNING: Skipping no-APA channel range: " << crn << endl;
322  continue;
323  }
324  std::istringstream ssapa(crn.substr(ipos,1));
325  ssapa >> iapa;
326  if ( iapa <= 0 ) {
327  // cout << myname << "ERROR: Unable to extract APA index from channel range " << crn << endl;
328  // continue;
329  // Allow no APA number for Iceberg.
330  iapa = -1;
331  }
332  }
333  const IndexRange& ran = pcrt->get(crn);
334  if ( ! ran.isValid() ) {
335  cout << myname << "WARNING: No channels taken from invalid channel range " << ran.name << endl;
336  cout << myname << "WARNING: Note that channel ranges are defined by tool channelRanges." << endl;
337  continue;
338  }
339  m_apacrns[iapa].push_back(crn);
340  for ( Index icha=ran.begin; icha<ran.end; ++icha ) {
341  bool keep = (keepChans.size() == 0 || keepChans.count(icha)) && skipChans.count(icha) == 0;
342  if ( keep ) m_apachsets[iapa].insert(icha);
343  }
344  if ( iapa == -1 ) {
345  nchaAll = ran.size();
346  break;
347  }
348  ++ncrn;
349  }
350 
351  // Display the APAs to be processed.
352  if ( m_apachsets.size() == 0 ) {
353  cout << myname << "WARNING: No channels will be processed." << endl;
354  } else if ( m_LogLevel >= 1 ) {
355  cout << myname << "The following APAs will be processed." << endl;
356  cout << myname << "-------------------------------------" << endl;
357  cout << myname << " APA # chan Channel ranges" << endl;
358  for ( ApaChannelSetMap::value_type iapa : m_apachsets ) {
359  cout << myname << setw(4);
360  if ( iapa.first == -1 ) cout << "all";
361  else cout << iapa.first;
362  cout << setw(9) << iapa.second.size() << " {";
363  bool first = true;
364  for ( Name crn : m_apacrns[iapa.first] ) {
365  if ( first ) first = false;
366  else cout << ", ";
367  cout << crn;
368  }
369  cout << "}" << endl;
370  }
371  cout << myname << "-------------------------------------" << endl;
372  }
373 
374  // Evaluate the maximum number of channels in the output containers.
375  AdcChannel nchaDefault = 2560;
376  AdcChannel nchaIn = 0;
377  AdcChannel nchaOut = 0;
378  if ( m_ApaChannelCounts.size() == 0 ) {
379  cout << myname << "WARNING: No APA channel counts have been provided. Using " << nchaDefault << endl;
380  } else {
381  nchaDefault = m_ApaChannelCounts.back();
382  }
383  for ( const ApaChannelSetMap::value_type& csmPair : m_apachsets ) {
384  int iapaSigned = csmPair.first;
385  AdcChannel nchaApa = nchaAll;
386  if ( iapaSigned >= 0 ) {
387  Index iapa = iapaSigned;
388  nchaApa = iapa < m_ApaChannelCounts.size() ? m_ApaChannelCounts[iapa] : nchaDefault;
389  }
390  nchaIn += nchaApa;
391  nchaOut += csmPair.second.size();
392  }
393  m_maxOutputDigitChannelCount = m_OutputDigitName.size() == 0 ? 0 : nchaIn;
394  m_maxOutputTimeStampChannelCount = m_OutputTimeStampName.size() == 0 ? 0 : nchaIn;
395  m_maxOutputWireChannelCount = m_OutputWireName.size() == 0 ? 0 : nchaOut;
396  if ( m_LogLevel >= 1 ) {
397  cout << myname << " Max # read channels: " << nchaIn << endl;
398  cout << myname << " Max # processed channels: " << nchaOut << endl;
399  cout << myname << " Max # output digits: " << m_maxOutputDigitChannelCount << endl;
400  cout << myname << " Max # output timestamps: " << m_maxOutputTimeStampChannelCount << endl;
401  cout << myname << " Max # output wires: " << m_maxOutputWireChannelCount << endl;
402  }
403 }
404 
405 //**********************************************************************
406 
408  const string myname = "DataPrepByApaModule::beginJob: ";
409  if ( m_LogLevel >= 2 ) cout << myname << "Starting job." << endl;
410  m_nproc = 0;
411  m_nskip = 0;
412 }
413 
414 //**********************************************************************
415 
417  const string myname = "DataPrepByApaModule::endJob: ";
418  if ( m_LogLevel >= 1 ) {
419  cout << myname << "# events processed: " << m_nproc << endl;
420  cout << myname << " # events skipped: " << m_nskip << endl;
421  }
422 }
423 
424 //**********************************************************************
425 
427  const string myname = "DataPrepByApaModule::produce: ";
428 
429  // Flag indicating that non-verbose info level messages should be logged.
430  bool logInfo = m_LogLevel >= 2;
431 
432  // Factor to convert event clock to ticks.
433  unsigned long triggerPerTick = 25;
434 
435  // Control flags.
436  bool skipAllEvents = false;
437  bool skipEventsWithCorruptDataDropped = false;
438 
439  // Decode the event time.
440  Timestamp beginTime = evt.time();
441  time_t itim = beginTime.timeHigh();
442  int itimrem = beginTime.timeLow();
443  // Older protoDUNE data has time in low field.
444  if ( itim == 0 && itimrem != 0 ) {
445  itimrem = itim;
446  itim = beginTime.timeLow();
447  }
448 
449  // Log event processing header.
450  if ( logInfo ) {
451  cout << myname << "Run " << evt.run();
452  if ( evt.subRun() ) cout << "-" << evt.subRun();
453  cout << ", event " << evt.event();
454  cout << ", nproc=" << m_nproc;
455  if ( m_nskip ) cout << ", nskip=" << m_nskip;
456  cout << endl;
457  if ( m_LogLevel >= 3 ) cout << myname << "Reading raw digits for producer, name: " << m_DigitProducer << ", " << m_DigitName << endl;
458  if ( evt.isRealData() ) {
459  TTimeStamp rtim(itim, itimrem);
460  string stim = string(rtim.AsString("s")) + " UTC";
461  cout << myname << "Real data event time: " << itim << " (" << stim << ")" << endl;
462  } else {
463  cout << myname << "Sim data event time: " << DuneTimeConverter::toString(beginTime) << endl;
464  }
465  }
466  if ( m_LogLevel >= 3 ) {
467  cout << myname << "Event time high, low: " << beginTime.timeHigh() << ", " << beginTime.timeLow() << endl;
468  }
469 
470  // Fetch the event trigger and timing clock.
471  AdcIndex trigFlag = 0;
472  AdcLongIndex timingClock = 0;
473  using TimeStampVector = std::vector<raw::RDTimeStamp>;
474  const TimeStampVector* ptims = nullptr;
475  if ( true ) {
476  art::InputTag itag1("timingrawdecoder", "daq");
477  auto htims = evt.getHandle<TimeStampVector>(itag1);
478  if ( htims ) {
479  ptims = &*htims;
480  if ( ptims->size() == 0 ) {
481  cout << myname << "No timing clocks found." << endl;
482  for ( unsigned int itim=0; itim<ptims->size() && itim<50; ++itim ) {
483  cout << myname << " " << ptims->at(itim).GetTimeStamp() << endl;
484  }
485  } else if ( ptims->size() > 1 ) {
486  cout << myname << "ERROR: Too many timing clocks: " << ptims->size() << endl;
487  for ( unsigned int itim=0; itim<ptims->size() && itim<50; ++itim ) {
488  cout << myname << " " << ptims->at(itim).GetTimeStamp() << endl;
489  }
490  } else {
491  const raw::RDTimeStamp& tim = ptims->at(0);
492  if ( logInfo ) cout << myname << "Timing clock: " << tim.GetTimeStamp() << endl;
493  timingClock = tim.GetTimeStamp();
494  // See https://twiki.cern.ch/twiki/bin/view/CENF/TimingSystemAdvancedOp#Reference_info
495  trigFlag = tim.GetFlags();
496  if ( m_LogLevel >= 2 ) cout << myname << "Trigger flag: " << trigFlag << " (";
497  bool isBeam = trigFlag == 0xc;
498  bool isCrt = trigFlag == 13;
499  bool isFake = trigFlag >= 0x8 && trigFlag <= 0xb;
500  if ( logInfo ) {
501  if ( isBeam ) cout << "Beam";
502  else if ( isCrt ) cout << "CRT";
503  else if ( isFake ) cout << "Fake";
504  else cout << "Unexpected";
505  cout << ")" << endl;
506  }
507  }
508  } else {
509  if ( logInfo ) {
510  cout << myname << "WARNING: Event timing clocks product not found." << endl;
511  }
512  }
513  }
514 
515  // Fetch beam information
516  float beamTof = 0.0; // Time of flight.
517  if ( m_BeamEventLabel.size() ) {
518  std::vector< art::Ptr<beam::ProtoDUNEBeamEvent> > beaminfo;
519  auto pdbeamHandle = evt.getHandle< std::vector<beam::ProtoDUNEBeamEvent> >(m_BeamEventLabel);
520  if ( pdbeamHandle ) {
521  art::fill_ptr_vector(beaminfo, pdbeamHandle);
522  if ( beaminfo.size() == 0 ) {
523  cout << myname << "WARNING: Beam event vector is empty." << endl;
524  } else {
525  if ( beaminfo.size() > 1 ) {
526  cout << myname << "WARNING: Beam event vector has size " << beaminfo.size() << endl;
527  }
528  AdcIndex beamTrigFlag = beaminfo[0]->GetTimingTrigger();
529  if ( beamTrigFlag != trigFlag ) {
530  if ( logInfo ) cout << myname << "Beam event and timing trigger flags differ: " << beamTrigFlag << " != " << trigFlag << endl;
531  } else if ( beamTrigFlag != 12 ) {
532  if ( logInfo ) cout << myname << "Beam event trigger is not beam: it is " << beamTrigFlag << endl;
533  //} else if ( ! beaminfo[0]->CheckIsMatched() ) {
534  // cout << myname << "Beam event is not matched." << endl;
535  } else if ( beaminfo[0]->GetTOFChan() == -1 ) {
536  if ( logInfo ) cout << myname << "Beam event does not have a TOF match." << endl;
537  } else {
538  int beamChan = beaminfo[0]->GetTOFChan();
539  beamTof = beaminfo[0]->GetTOF();
540  if ( logInfo ) cout << myname << "Beam event TOF[" << beamChan << "]: " << beamTof << endl;
541  }
542  const std::vector<double>& momenta = beaminfo[0]->GetRecoBeamMomenta();
543  if ( logInfo ) {
544  cout << myname << "Beam momenta:";
545  for ( double pval : momenta ) cout << " " << pval;
546  cout << endl;
547  }
548  }
549  } else {
550  if ( logInfo ) cout << myname << "Beam event data product not found: " << m_BeamEventLabel << endl;
551  }
552  }
553 
554  // Create containers for data to be written to the event data store.
555  using DigitVector = std::vector<raw::RawDigit>;
556  std::unique_ptr<DigitVector> pdigitsAll;
558  pdigitsAll.reset(new DigitVector);
559  pdigitsAll->reserve(m_maxOutputDigitChannelCount);
560  }
561  std::unique_ptr<TimeStampVector> ptimsAll;
563  ptimsAll.reset(new TimeStampVector);
564  ptimsAll->reserve(m_maxOutputTimeStampChannelCount);
565  }
566  using WireVector = std::vector<recob::Wire>;
567  std::unique_ptr<std::vector<recob::Wire>> pwires;
569  pwires.reset(new WireVector);
570  pwires->reserve(m_maxOutputWireChannelCount);
571  }
572  using StatVector = std::vector<raw::RDStatus>;
573  std::unique_ptr<StatVector> pstatusAll;
575  pstatusAll.reset(new StatVector);
576  pstatusAll->reserve(m_maxOutputDigitChannelCount);
577  }
578 
579  // Notify data preparation service of start of event.
580  DuneEventInfo devt;
581  devt.run = evt.run();
582  devt.subRun = evt.subRun();
583  devt.event = evt.event();
584  devt.triggerClock = timingClock;
585  devt.triggerTick0 = timingClock/triggerPerTick;
586  devt.time = itim;
587  devt.timerem = itimrem;
588  int bstat = m_pRawDigitPrepService->beginEvent(devt);
589  if ( bstat ) cout << myname << "WARNING: Event initialization failed." << endl;
591 
592  // Loop over channel ranges.
594  if ( ptm == nullptr ) {
595  cout << myname << "ERROR: Tool manager not found." << endl;
596  return;
597  }
598  const IndexRangeTool* pcrt = ptm->getShared<IndexRangeTool>("channelRanges");
599  if ( pcrt == nullptr ) {
600  cout << myname << "ERROR: Channel range tool channelRanges not found." << endl;
601  return;
602  }
603  // Loop over APAs or all.
604  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
605  for ( const ApaChannelSetMap::value_type& csmPair : m_apachsets ) {
606  int iapa = csmPair.first;
607  ostringstream ssapa;
608  if ( iapa == -1 ) {
609  ssapa << "All";
610  } else {
611  ssapa << "APA " << iapa;
612  }
613  Name sapa = ssapa.str();
614  std::vector<int> apas = {iapa};
615  // Fetch the digits, clocks and status for the channel range.
616  if ( logInfo ) cout << myname << "Fetching digits and clocks for " << sapa << "." << endl;
617  using StatVector = std::vector<raw::RDStatus>;
618  TimeStampVector timsCrn;
619  StatVector statsCrn;
620  DigitVector digitsCrn;
621  // Fetch the data for this APA.
622  int decodeStat = m_pDecoderTool->
623  retrieveDataForSpecifiedAPAs(evt, digitsCrn, timsCrn, statsCrn, apas);
624  if ( m_LogLevel >= 3 ) { // Decoder tool can return any value for success
625  cout << myname << "INFO: Decoder tool for APA " << iapa << " returned " << decodeStat << endl;
626  }
627  if ( logInfo ) {
628  cout << myname << " " << sapa << " digit count from tool: " << digitsCrn.size() << endl;
629  cout << myname << " " << sapa << " stats count from tool: " << statsCrn.size() << endl;
630  cout << myname << " " << sapa << " clock count from tool: " << timsCrn.size() << endl;
631  }
632  // Check read status.
633  string srdstat;
634  bool skipEvent = skipAllEvents;
635  if ( statsCrn.size() != 1 ) cout << myname << "WARNING: Read status has unexpected size: "
636  << statsCrn.size() << endl;
637  const RDStatus rdstat = statsCrn.at(0);
638  if ( false ) {
639  cout << myname << "Raw data status: " << rdstat.GetStatWord();
640  if ( rdstat.GetCorruptDataDroppedFlag() ) cout << " (Corrupt data was dropped.)";
641  if ( rdstat.GetCorruptDataKeptFlag() ) cout << " (Corrupt data was retained.)";
642  cout << endl;
643  }
644  cout << myname << "Raw data read status: " << std::to_string(rdstat.GetStatWord()) << endl;
645  srdstat = "rdstat=" + std::to_string(rdstat.GetStatWord());
646  skipEvent |= skipEventsWithCorruptDataDropped && rdstat.GetCorruptDataDroppedFlag();
647  // Check the channel clocks from the tool.
648  vector<AdcLongIndex> channelClocks;
649  vector<ULong64_t> tzeroClockCandidates; // Candidates for t0 = 0.
650  float trigTickOffset = -500.5;
651  using ClockCounter = std::map<ULong64_t, AdcIndex>;
652  using ClockDiffs = std::map<ULong64_t, long>;
653  using ClockTickDiffs = std::map<ULong64_t, float>;
654  using ClockMessages = std::map<ULong64_t, Name>;
655  using NtickCounter = std::map<AdcIndex, AdcIndex>;
656  ClockCounter clockCounts;
657  ClockDiffs clockDiffs;
658  ClockTickDiffs tickDiffs;
659  ClockMessages clockMessages;
660  NtickCounter ntickCounter;
661  ULong64_t chClock = 0;
662  ULong64_t maxdiff = 99999999; // 2 sec
663  float tickdiff = maxdiff;
664  AdcIndex nskipEmpty = 0;
665  unsigned int ntim = timsCrn.size();
666  unsigned int ndigi = digitsCrn.size();
667  if ( ntim > 0 && ntim != ndigi ) {
668  cout << "ERROR: Channel clock count differs from digit count: " << ntim << " !' " << ndigi << endl;
669  abort();
670  }
671  for ( unsigned int idig=0; idig<ntim; ++idig ) {
672  raw::RDTimeStamp chts = timsCrn[idig];
673  raw::RawDigit& dig = digitsCrn[idig];
674  AdcIndex ntick = dig.Samples();
675  if ( ntickCounter.count(ntick) == 0 ) ntickCounter[ntick] = 1;
676  else ++ntickCounter[ntick];
677  if ( ntick == 0 ) { // Do not check clocks for empty digits
678  ++nskipEmpty;
679  continue;
680  }
681  chClock = chts.GetTimeStamp();
682  channelClocks.push_back(chClock);
683  bool sign = chClock > timingClock;
684  ULong64_t chClockAbsDiff = sign ? chClock - timingClock
685  : timingClock - chClock;
686  bool badDiff = chClockAbsDiff > maxdiff;
687  if ( badDiff ) chClockAbsDiff = maxdiff;
688  long chClockDiff = sign ? chClockAbsDiff : -chClockAbsDiff;
689  tickdiff = chClockDiff/triggerPerTick;
690  bool nearTrigger = fabs(tickdiff - trigTickOffset) < 1.0;
691  if ( clockCounts.find(chClock) == clockCounts.end() ) {
692  clockCounts[chClock] = 1;
693  clockDiffs[chClock] = chClockDiff;
694  tickDiffs[chClock] = tickdiff;
695  Name msg;
696  if ( chClock == 0 ) msg = "Channel clock is zero.";
697  else if ( badDiff ) msg = "Channel clock is very far from timing clock.";
698  clockMessages[chClock] = msg;
699  if ( nearTrigger ) tzeroClockCandidates.push_back(chClock);
700  } else {
701  ++clockCounts[chClock];
702  }
703  if ( ! nearTrigger && timingClock != 0 ) {
704  if ( m_LogLevel >= 3 ) {
705  cout << myname << "WARNING: Channel timing difference: " << chClockDiff
706  << " (" << tickdiff << " ticks)." << endl;
707  }
708  }
709  }
710  if ( clockCounts.size() > 1 ) {
711  if ( logInfo ) {
712  cout << myname << "WARNING: Channel clocks for " << sapa << " are not consistent." << endl;
713  cout << myname << "WARNING: Clock ticks count" << endl;
714  for ( ClockCounter::value_type iclk : clockCounts ) {
715  ULong64_t chClock = iclk.first;
716  AdcIndex count = iclk.second;
717  long chClockDiff = clockDiffs[chClock];
718  float tickdiff = tickDiffs[chClock];
719  Name msg = clockMessages[chClock];
720  if ( logInfo ) {
721  cout << myname << "WARNING:" << setw(10) << chClockDiff << setw(10) << tickdiff
722  << setw(8) << count;
723  if ( msg.size() ) cout << " " << msg;
724  cout << endl;
725  }
726  }
727  if ( nskipEmpty ) {
728  cout << myname << "WARNING: No ADC samples:" << setw(8) << nskipEmpty << endl;
729  }
730  }
731  } else if ( clockCounts.size() == 1 ) {
732  if ( logInfo ) {
733  cout << myname << "Channel clocks for " << sapa << " are consistent with an offset of "
734  << tickdiff << " ticks";
735  if ( nskipEmpty ) cout << " (" << nskipEmpty << " channels skipped)";
736  cout << "." << endl;
737  }
738  } else if ( ntim > 0 ) {
739  // This should never happen.
740  if ( logInfo ) cout << myname << "WARNING: Channel clocks not found." << endl;
741  }
742  // Find the primary (most common) tick count.
743  AdcIndex ntickMaxCount = 0; // tick count with the most channels
744  AdcIndex nchanMaxCount = 0; // # channels with the tick count
745  for ( NtickCounter::value_type ent : ntickCounter ) {
746  if ( ent.second > nchanMaxCount ) {
747  ntickMaxCount = ent.first;
748  nchanMaxCount = ent.second;
749  }
750  }
751  AdcIndex ntickPrimary = ntickMaxCount;
752  // Evaluate the allowed range of tick counts.
753  AdcIndex minTickCount = 0;
754  AdcIndex maxTickCount = 0;
755  if ( m_DeltaTickCount == 0 ) {
756  minTickCount = ntickPrimary;
757  maxTickCount = ntickPrimary;
758  } else if ( m_DeltaTickCount > 0 ) {
759  float tcmin = ntickPrimary*(1.0 - m_DeltaTickCount);
760  float tcmax = ntickPrimary*(1.0 + m_DeltaTickCount);
761  minTickCount = tcmin > 0.0 ? tcmin : 0;
762  maxTickCount = tcmax > 0.0 ? tcmax : 0;
763  cout << myname << "Allowed tick count range is (" << minTickCount << ", " << maxTickCount << ")" << endl;
764  }
765  // Remove out-of-range tick counts.
766  // Those channels will not be processed.
767  if ( maxTickCount > 0 ) {
768  NtickCounter::iterator ient = ntickCounter.begin();
769  while ( ient != ntickCounter.end() ) {
770  AdcIndex ntick = ient->first;
771  if ( ntick < minTickCount || ntick > maxTickCount ) {
772  ient = ntickCounter.erase(ient);
773  } else {
774  ++ient;
775  }
776  }
777  }
778  // Display tick counts.
779  if ( logInfo && ntickCounter.size() ) {
780  if ( ntickCounter.size() == 1 ) {
781  cout << myname << "Tick count for all channels is " << ntickCounter.begin()->first << endl;
782  } else {
783  //string slev = "WARNING: ";
784  string slev = "";
785  cout << myname << slev << "Retaining inconsistent tick counts:" << endl;
786  cout << myname << slev << " Ntick Nchan" << endl;
787  for ( NtickCounter::value_type ent : ntickCounter ) {
788  cout << myname << slev << setw(8) << ent.first << setw(8) << ent.second << endl;
789  }
790  }
791  }
792  // Build the AdcChannelData objects.
793  AdcChannelDataMap datamap;
794  if ( m_LogLevel >= 3 ) {
795  cout << myname << "# digits read for " << sapa << ": " << digitsCrn.size() << endl;
796  }
797  // Create the transient data map and copy the digits there.
798  const ApaChannelSet& csmKeepChans = csmPair.second;
799  unsigned int nproc = 0;
800  unsigned int nkeep = 0;
801  unsigned int nskip = 0;
802  unsigned int nempty = 0;
803  unsigned int nbadtc = 0;
804  for ( unsigned int idig=0; idig<ndigi; ++idig ) {
805  raw::RawDigit& dig = digitsCrn[idig];
806  AdcChannel chan = dig.Channel();
807  if ( csmKeepChans.count(chan) == 0 ) {
808  ++nskip;
809  continue;
810  }
811  // 02jan2020: Skip empty channels (Redmine 23811).
812  AdcIndex ntick = dig.Samples();
813  if ( ntick == 0 ) {
814  ++nempty;
815  if ( m_SkipEmptyChannels ) {
816  ++nskip;
817  continue;
818  }
819  }
820  // 06jan202: Skip channels with too few or too many ticks.
821  if ( maxTickCount > minTickCount ) {
822  if ( ntick < minTickCount || ntick > maxTickCount ) {
823  ++nbadtc;
824  ++nskip;
825  continue;
826  }
827  }
828  // Create AdcChannelData for this channel.
829  auto its = datamap.emplace(chan, AdcChannelData());
830  if ( ! its.second ) {
831  cout << myname << "ERROR: Ignoring duplicate digit for channel " << chan << "." << endl;
832  ++nskip;
833  continue;
834  }
835  AdcChannelData& acd = its.first->second;
836  ++nproc;
837  // Fetch the channel status.
838  Index chanStat = AdcChannelStatusGood;
839  if ( m_pChannelStatusProvider != nullptr ) {
840  if ( m_pChannelStatusProvider->IsNoisy(chan) ) chanStat = AdcChannelStatusNoisy;
841  if ( m_pChannelStatusProvider->IsBad(chan) ) chanStat = AdcChannelStatusBad;
842  }
843  // Fetch the online ID.
845  AdcChannel fembChannel = AdcChannelData::badIndex();
846  if ( m_onlineChannelMapTool ) {
847  unsigned int ichOn = m_onlineChannelMapTool->get(chan);
848  if ( ichOn != IndexMapTool::badIndex() ) {
849  fembID = ichOn/128;
850  fembChannel = ichOn % 128;
851  }
852  }
853  // Build the channel data.
854  acd.setEventInfo(pevt);
855  acd.setChannelInfo(chan, fembID, fembChannel, chanStat);
856  acd.digitIndex = idig;
857  acd.digit = &dig;
858  if ( channelClocks.size() > idig ) {
859  acd.channelClock = channelClocks[idig];
860  }
861  acd.metadata["ndigi"] = ndigi;
862  if ( m_BeamEventLabel.size() ) {
863  acd.metadata["beamTof"] = beamTof;
864  }
865  ++nkeep;
866  }
867  if ( logInfo ) {
868  cout << myname << " " << sapa << " # input digits: " << ndigi << endl;
869  cout << myname << " " << sapa << " # channels selected: " << nkeep;
870  if ( ! m_SkipEmptyChannels && nempty ) cout << " (" << nempty << " empty)";
871  cout << endl;
872  cout << myname << " " << sapa << " # channels skipped: " << nskip;
873  bool havebad = false;
874  if ( m_SkipEmptyChannels && nempty ) {
875  cout << " (" << nempty << " empty";
876  havebad = true;
877  }
878  if ( nbadtc ) {
879  cout << (havebad ? ", " : " (") << nbadtc << " bad tick count";
880  havebad = true;
881  }
882  if ( havebad ) cout << ")";
883  cout << endl;
884  cout << myname << " " << sapa << " # channels to be processed: " << nproc << endl;
885  }
886 
887  Index nacd = datamap.size();
888  if ( nacd == 0 ) {
889  if ( m_LogLevel >= 3 ) {
890  if ( logInfo ) cout << myname << "Skipping empty data map." << endl;
891  }
892  continue;
893  }
894 
895  // Make sure created wires are deleted or stored in event.
896  WireVector wiresTmp;
897  WireVector* pwiresCrn = pwires ? pwires.get() : &wiresTmp;
898 
899  // Use the data preparation service to build the wires.
900  if ( m_LogLevel >= 3 ) {
901  cout << myname << "Preparing " << nacd << " channel" << (nacd == 1 ? "" : "s") << "." << endl;
902  }
903  int rstat = m_pRawDigitPrepService->prepare(clockData, datamap, pwiresCrn, nullptr);
904  if ( rstat != 0 ) {
905  cout << myname << "ERROR: Data preparation service returned error " << rstat;
906  continue;
907  }
908 
909  // Transfer the larsoft digits to the output container.
910  if ( pdigitsAll ) for ( raw::RawDigit& dig : digitsCrn ) pdigitsAll->emplace_back(dig);
911  if ( ptimsAll ) for ( raw::RDTimeStamp& tst : timsCrn ) ptimsAll->emplace_back(tst);
912  if ( pstatusAll ) for ( raw::RDStatus& stat : statsCrn ) pstatusAll->emplace_back(stat);
913 
914  } // End loop over channel ranges.
915 
916  // Notify data preparation service of end of event.
917  int estat = m_pRawDigitPrepService->endEvent(devt);
918  if ( estat ) cout << myname << "WARNING: Event finalization failed." << endl;
919 
920  // Record wires and associations in the event.
921  if ( pwires ) {
922  if ( logInfo ) cout << myname << "Created wire count: " << pwires->size() << endl;
923  if ( pwires->size() == 0 ) cout << myname << "WARNING: No wires made for this event." << endl;
924  evt.put(std::move(pwires), m_OutputWireName);
925  } else {
926  if ( logInfo ) cout << myname << "Wire output was not requested." << endl;
927  }
928 
929  // Record decoder containers.
930  if ( m_OutputDigitName.size() ) {
931  if ( pdigitsAll ) {
932  if ( logInfo ) cout << myname << "Created digit count: " << pdigitsAll->size() << endl;
933  evt.put(std::move(pdigitsAll), m_OutputDigitName);
934  } else {
935  cout << myname << "WARNING: Digits are not written because container was not created." << endl;
936  }
937  } else {
938  if ( logInfo ) cout << myname << "Digit output was not requested." << endl;
939  }
940 
941  if ( m_OutputTimeStampName.size() ) {
942  if ( ptimsAll ) {
943  if ( logInfo ) cout << myname << "Created time stamp count: " << ptimsAll->size() << endl;
944  evt.put(std::move(ptimsAll), m_OutputTimeStampName);
945  } else {
946  cout << myname << "WARNING: Output time stamp container was not created." << endl;
947  }
948  } else {
949  if ( logInfo ) cout << myname << "Time stamp output was not requested." << endl;
950  }
951 
952  if ( m_OutputDigitName.size() ) {
953  if ( pstatusAll ) {
954  if ( logInfo ) cout << myname << "Created digit count: " << pstatusAll->size() << endl;
955  evt.put(std::move(pstatusAll), m_OutputDigitName);
956  } else {
957  cout << myname << "WARNING: Output status container was not created." << endl;
958  }
959  } else {
960  if ( logInfo ) cout << myname << "Status output was not requested." << endl;
961  }
962 
963  ++m_nproc;
964  return;
965 }
virtual bool IsBad(raw::ChannelID_t channel) const =0
Returns whether the specified channel is bad in the current run.
void reconfigure(fhicl::ParameterSet const &p)
void produce(art::Event &evt)
intermediate_table::iterator iterator
LongIndex triggerTick0
Definition: DuneEventInfo.h:28
constexpr std::uint32_t timeLow() const
Definition: Timestamp.h:29
EventNumber_t event() const
Definition: DataViewImpl.cc:85
ULong64_t Samples() const
Number of samples in the uncompressed ADC data.
Definition: RawDigit.h:213
uint16_t GetFlags() const
Definition: RDTimeStamp.h:46
Collection of charge vs time digitized from a single readout channel.
Definition: RawDigit.h:69
void msg(const char *fmt,...)
Definition: message.cpp:107
const AdcIndex AdcChannelStatusNoisy
Definition: AdcTypes.h:48
static std::string toString(art::Timestamp tart)
Handle< PROD > getHandle(SelectorBase const &) const
Definition: DataViewImpl.h:382
std::string string
Definition: nybbler.cc:12
Index begin
Definition: IndexRange.h:34
virtual bool IsNoisy(raw::ChannelID_t channel) const =0
Returns whether the specified channel is noisy in the current run.
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.h:20
constexpr std::uint32_t timeHigh() const
Definition: Timestamp.h:34
ChannelID_t Channel() const
DAQ channel this raw data was read from.
Definition: RawDigit.h:212
struct vector vector
Index size() const
Definition: IndexRange.h:88
bool isValid() const
Definition: IndexRange.h:94
std::shared_ptr< const EventInfo > EventInfoPtr
Index end
Definition: IndexRange.h:35
const raw::RawDigit * digit
virtual int prepare(detinfo::DetectorClocksData const &clockData, AdcChannelDataMap &prepdigs, std::vector< recob::Wire > *pwires=nullptr, WiredAdcChannelDataMap *pwiredData=nullptr) const =0
std::vector< AdcChannel > AdcChannelVector
Definition: AdcTypes.h:51
Name name
Definition: IndexRange.h:32
std::vector< Name > NameVector
std::unique_ptr< IndexMapTool > m_onlineChannelMapTool
bool GetCorruptDataDroppedFlag() const
Definition: RDStatus.h:31
std::unique_ptr< PDSPTPCDataInterfaceParent > m_pDecoderTool
bool isRealData() const
void setChannelInfo(ChannelInfoPtr pchi)
std::map< int, ApaChannelSet > ApaChannelSetMap
Timestamp time() const
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
const lariov::ChannelStatusProvider * m_pChannelStatusProvider
std::map< int, NameVector > ApaNamesMap
def move(depos, offset)
Definition: depos.py:107
T get(std::string const &key) const
Definition: ParameterSet.h:271
void setEventInfo(EventInfoPtr pevi)
RawDigitPrepService * m_pRawDigitPrepService
p
Definition: test.py:223
SubRunNumber_t subRun() const
Definition: DataViewImpl.cc:78
ProductID put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: DataViewImpl.h:686
Class providing information about the quality of channels.
RunNumber_t run() const
Definition: DataViewImpl.cc:71
virtual IndexRangeGroup get(Name nam) const =0
unsigned int AdcIndex
Definition: AdcTypes.h:15
Name m_OutputWireName
Second field in full label for the output wire container.
ULong64_t GetTimeStamp() const
Definition: RDTimeStamp.h:42
Q_EXPORT QTSManip setw(int w)
Definition: qtextstream.h:331
virtual int beginEvent(const DuneEventInfo &) const
int sign(double val)
Definition: UtilFunc.cxx:104
pval
Definition: tracks.py:168
unsigned long AdcLongIndex
Definition: AdcTypes.h:16
RangeVector ranges
std::vector< art::Ptr< recob::Wire > > WireVector
unsigned int AdcChannel
Definition: AdcTypes.h:50
static Index badIndex()
Definition: IndexMapTool.h:18
std::optional< T > get_if_present(std::string const &key) const
Definition: ParameterSet.h:224
const AdcIndex AdcChannelStatusGood
Definition: AdcTypes.h:46
virtual int endEvent(const DuneEventInfo &) const
AdcLongIndex channelClock
Interface for experiment-specific channel quality info provider.
bool GetCorruptDataKeptFlag() const
Definition: RDStatus.h:32
Declaration of basic channel signal object.
std::map< AdcChannel, AdcChannelData > AdcChannelDataMap
const AdcIndex AdcChannelStatusBad
Definition: AdcTypes.h:47
AdcIndex digitIndex
TCEvent evt
Definition: DataStructs.cxx:7
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
Interface for experiment-specific service for channel quality info.
static Index badIndex()
bool isValid() const
unsigned int GetStatWord() const
Definition: RDStatus.h:33
std::string to_string(ModuleType const mt)
Definition: ModuleType.h:34
DataPrepByApaModule(fhicl::ParameterSet const &pset)
static DuneToolManager * instance(std::string fclname="", int dbg=1)
T * getShared(std::string name)
QTextStream & endl(QTextStream &s)
LongIndex triggerClock
Definition: DuneEventInfo.h:27