Dune35tNoiseRemovalService_service.cc
Go to the documentation of this file.
1 // Dune35tNoiseRemovalService_service.cc
2 
4 #include <iostream>
5 #include <sstream>
6 #include <iomanip>
7 #include <algorithm>
10 #include "lbne-raw-data/Services/ChannelMap/ChannelMapService.h"
11 
12 using std::vector;
13 using std::string;
14 using std::ostream;
15 using std::cout;
16 using std::endl;
17 using std::ostringstream;
18 using std::setw;
19 using art::ServiceHandle;
20 
21 //**********************************************************************
22 
25 : m_LogLevel(1), m_nwarn(0) {
26  const string myname = "Dune35tNoiseRemovalService::ctor: ";
27  pset.get_if_present<int>("LogLevel", m_LogLevel);
28  m_GroupingFlag = pset.get<int>("GroupingFlag");
29  m_SkipStuckCodes = pset.get<bool>("SkipStuckCodes");
30  m_SkipSignals = pset.get<bool>("SkipSignals");
31  m_CorrectStuckCodes = pset.get<bool>("CorrectStuckCodes");
32  m_ShowGroups = pset.get<int>("ShowGroups");
33  m_ShowGroupsChannel = pset.get<int>("ShowGroupsChannel");
34  // Get services.
37  // Build the map of offline channel labels.
38  const unsigned int nchan = m_pGeometry->Nchannels();
39  for ( unsigned int chanon=0; chanon<nchan; ++chanon ) {
40  unsigned int chanoff = m_pChannelMap->Offline(chanon);
41  geo::View_t gview = m_pGeometry->View(chanoff);
42  unsigned int iview = gview;
43 #if 0
44  // Switch to this when larcore GeometryCore supports ROPs.
45  // https://cdcvs.fnal.gov/redmine/issues/9264
46  readout::ROPID ropid = m_pGeometry->ChannelToROP(chanoff);
47  unsigned int rop = ropid.ROP;
48  unsigned int chan0 = m_pGeometry->FirstChannelInROP(ropid);
49  unsigned int ropchan = chanoff - chan0;
50 #else
51  // Hard wire 35 ton convention.
52  unsigned int ropchan = 999;
53  {
54  geo::WireID wid = m_pGeometry->ChannelToWire(chanoff).at(0);
55  unsigned int tpcoff = wid.TPC;
56  if ( gview == geo::kW && tpcoff%2 ) iview = 3;
57  unsigned int apachan = chanoff%512;
58  if ( iview == 0 ) ropchan = apachan;
59  if ( iview == 1 ) ropchan = apachan - 144;
60  if ( iview == 2 ) ropchan = apachan - 288;
61  if ( iview == 3 ) ropchan = apachan - 400;
62  }
63 #endif
64  static string sview[4] = {"u", "v", "z1", "z2"};
65  unsigned int apa = m_pChannelMap->APAFromOnlineChannel(chanon);
66  ostringstream ssrop;
67  ssrop << apa << sview[iview];
68  string srop = ssrop.str(); // Name for the ROP
69  ssrop << "-";
70  if ( ropchan < 100 ) ssrop << "0";
71  if ( ropchan < 10 ) ssrop << "0";
72  ssrop << ropchan;
73  string sropchan = ssrop.str(); // Name for the ROP channel
74  m_sRopChannelMap[chanon] = sropchan;
75  }
76  // Build the channel groupings.
77  if ( m_GroupingFlag == 0 ) { //group by regulators
78  m_GroupChannels.resize(3);
79  for (size_t i = 0; i<3; ++i){
80  m_GroupChannels[i].resize(32);
81  }
82  for (unsigned int i=0; i<nchan; ++i){//online channels
83  unsigned int plane = m_pChannelMap->PlaneFromOnlineChannel(i);
84  unsigned int rce = m_pChannelMap->RCEFromOnlineChannel(i);
85  unsigned int regulator = m_pChannelMap->RegulatorFromOnlineChannel(i);
86  m_GroupChannels[plane][rce*2+regulator].push_back(i);
87  }
88  } else if ( m_GroupingFlag == 1 ) {
89  m_GroupChannels.resize(3);
90  for (size_t i = 0; i<3; ++i){
91  m_GroupChannels[i].resize(128);
92  }
93  if ( m_LogLevel >= 2 ) {
94  cout << myname << "Online channel: plane:RCE:ASIC ROP-ropchan" << endl;
95  }
96  for ( unsigned int chanon=0; chanon<nchan; ++chanon ) { //online channels
97  unsigned int plane = m_pChannelMap->PlaneFromOnlineChannel(chanon);
98  unsigned int rce = m_pChannelMap->RCEFromOnlineChannel(chanon);
99  unsigned int asic = m_pChannelMap->ASICFromOnlineChannel(chanon);
100  if ( m_LogLevel >= 2 ) {
101  cout << myname << chanon << ": " << plane << ":" << rce << ":" << asic
102  << " " << m_sRopChannelMap.find(chanon)->second
103  << endl;
104  }
105  m_GroupChannels[plane][rce*8+asic].push_back(chanon);
106  }
107  } else {
108  cout << myname << "ERROR: Invalid GroupingFlag: " << m_GroupingFlag << endl;
109  return;
110  }
111  print(cout, myname);
112 }
113 
114 //**********************************************************************
115 
117  const string myname = "Dune35tNoiseRemovalService:update: ";
118  if ( datamap.size() == 0 ) {
119  cout << myname << "WARNING: No channels found." << endl;
120  return 1;
121  }
122  unsigned int nsig = 0;
123  bool first = true;
124  for ( const AdcChannelDataMap::value_type& ent : datamap ) {
125  const AdcChannelData& data = ent.second;
126  if ( first ) nsig = data.samples.size();
127  if ( data.samples.size() != nsig ) {
128  cout << myname << "WARNING: Channels have inconsistent sample counts." << endl;
129  return 2;
130  }
131  }
132  if ( nsig == 0 ) {
133  cout << myname << "WARNING: No ADC samples found." << endl;
134  return 3;
135  }
136  if ( m_GroupChannels.size() == 0 ) {
137  cout << myname << "ERROR: Channel groupings have not been defined." << endl;
138  return 4;
139  }
140  if ( m_LogLevel >= 3 ) {
141  cout << myname << "Entering..." << endl;
142  cout << myname << " # Channels: " << datamap.size() << endl;
143  cout << myname << " # Ticks: " << nsig << endl;
144  }
145  // loop through time slices
146  for ( unsigned int isig=0; isig<nsig; ++isig) {
147  // Loop over wire orientations.
148  for ( size_t iori=0; iori<m_GroupChannels.size(); ++iori ){
149  // Loop over channel groups.
150  for ( size_t igrp=0; igrp<m_GroupChannels[iori].size(); ++igrp ){
151  const AdcChannelVector& chans = m_GroupChannels[iori][igrp];
152  if ( m_LogLevel >= 3 && isig==0 ) {
153  cout << myname << setw(3) << iori << setw(4) << igrp << " " << chans.size() << endl;
154  }
155  // Loop over the channels for this time slice, orientation and group
156  vector<double> corrVals;
157  for ( AdcChannel chan : chans ) {
158  unsigned int offlineChan = m_pChannelMap->Offline(chan);
159  if ( m_LogLevel >= 4 && isig==0 ) {
160  cout << myname << "Channel on/off: " << chan << "/" << offlineChan << endl;
161  }
162  AdcChannelDataMap::const_iterator iacd = datamap.find(offlineChan);
163  if ( iacd == datamap.end() ) continue;
164  const AdcChannelData& acd = iacd->second;
165  AdcSignal sig = acd.samples[isig];
166  AdcFlag flag = acd.flags.size() ? acd.flags[isig] : AdcGood;
167  // Skip stuck codes.
168  // The original code also skipped samples where the raw count or pedestal was less than 10.
169  if ( m_SkipStuckCodes && ( flag==AdcStuckOff || flag==AdcStuckOn ) ) continue;
170  // Skip signals.
171  bool isSignal = false;
172  if ( m_SkipSignals ) {
173  if ( acd.signal.size() > isig ) {
174  isSignal = acd.signal[isig];
175  } else if ( m_nwarn < 100 ) {
176  ++m_nwarn;
177  cout << myname << "WARNING: Signal skip requested without signal finding." << endl;
178  }
179  }
180  if ( m_SkipSignals && isSignal ) continue;
181  // Add current sample to correction.
182  corrVals.push_back(sig);
183  }
184  // Evaluate the correction as the median of the selected signals.
185  unsigned int corrValSize = corrVals.size();
186  sort(corrVals.begin(),corrVals.end());
187  double correction = 0;
188  if ( corrValSize < 2 )
189  correction = 0.0;
190  else if ( (corrValSize % 2) == 0 )
191  correction = (corrVals[corrValSize/2] + corrVals[(corrValSize/2)-1])/2.0;
192  else
193  correction = corrVals[(corrValSize-1)/2];
194  // Loop over the channels for this group
195  for ( AdcChannel chan : chans ) {
196  unsigned int offlineChan = m_pChannelMap->Offline(chan);
197  AdcChannelDataMap::iterator iacd = datamap.find(offlineChan);
198  if ( iacd == datamap.end() ) continue;
199  AdcChannelData& acd = iacd->second;
200  AdcSignal& sig = acd.samples[isig];
201  AdcFlag flag = acd.flags.size() ? acd.flags[isig] : AdcGood;
202  if ( m_SkipStuckCodes && ( flag==AdcStuckOff || flag==AdcStuckOn ) ) continue;
203  if ( datamap.count(offlineChan) == 0 ) continue;
204  if ( !m_CorrectStuckCodes && ( flag==AdcStuckOff || flag==AdcStuckOn ) ) continue;
205  sig -= correction;
206  } // end loop over channels in the group
207  } // end loop over groups
208  } // end loop over orientations
209  } // end loop over time samples
210  return 0;
211 }
212 
213 //**********************************************************************
214 
216 print(ostream& out, string prefix) const {
217  out << prefix << "Dune35tNoiseRemovalService:" << endl;
218  out << prefix << " LogLevel: " << m_LogLevel << endl;
219  out << prefix << " GroupingFlag: " << m_GroupingFlag << endl;
220  out << prefix << " SkipStuckCodes: " << m_SkipStuckCodes << endl;
221  out << prefix << " CorrectStuckCodes: " << m_CorrectStuckCodes << endl;
222  out << prefix << " ShowGroups: " << m_ShowGroups << endl;
223  out << prefix << " ShowGroupsChannel: " << m_ShowGroupsChannel << endl;
224  if ( m_ShowGroups ) {
225  int wori = 10;
226  int wgrp = 10;
227  int wcha = 5;
228  if ( m_ShowGroups == 1 ) {
229  out << prefix << setw(wori) << "Orient." << setw(wgrp) << "Group" << ": " << " Online channels" << endl;
230  } else if ( m_ShowGroups == 2 ) {
231  out << prefix << setw(wgrp+4) << "Group" << ": ";
232  if ( m_ShowGroupsChannel == 0 ) out << " Channels" << endl;
233  else if ( m_ShowGroupsChannel == 1 ) out << " Online channels" << endl;
234  else if ( m_ShowGroupsChannel <= 4 ) out << " Offline channels" << endl;
235  }
236  vector<vector<string>> gochans; // Store channel lists as strings ordering group before orient
237  unsigned int nori = m_GroupChannels.size();
238  for ( unsigned int iori=0; iori<nori; ++iori ) {
239  const std::vector<AdcChannelVector>& gchans = m_GroupChannels[iori];
240  for ( unsigned int igrp=0; igrp<gchans.size(); ++igrp ) {
241  const AdcChannelVector& chans = gchans[igrp];
242  if ( chans.size() == 0 ) continue;
243  std::vector<string> schans;
244  bool usesetw = false;
245  bool addspace = false;
246  for ( AdcChannel chan : chans ) {
247  if ( m_ShowGroupsChannel == 0 ) {
248  schans.push_back(".");
249  } else if ( m_ShowGroupsChannel == 1 ) {
250  ostringstream sout;
251  sout << chan;
252  schans.push_back(sout.str());
253  usesetw = true;
254  } else if ( m_ShowGroupsChannel == 2 ) {
255  ostringstream sout;
256  sout << m_pChannelMap->Offline(chan);
257  schans.push_back(sout.str());
258  usesetw = true;
259  } else if ( m_ShowGroupsChannel == 3 ) {
260  schans.push_back(m_sRopChannelMap.find(chan)->second);
261  addspace = true;
262  } else if ( m_ShowGroupsChannel == 4 ) {
263  string schan = m_sRopChannelMap.find(chan)->second;
264  string::size_type pos = schan.find("z1-");
265  if ( pos != string::npos ) schan.replace(pos, 2, "z");
266  pos = schan.find("z2-");
267  if ( pos != string::npos ) schan.replace(pos, 2, "Z");
268  schans.push_back(schan);
269  addspace = true;
270  }
271  }
272  std::sort(schans.begin(), schans.end());
273  ostringstream sschanlist;
274  for ( string schan : schans ) {
275  if ( usesetw ) sschanlist << setw(wcha);
276  else if ( addspace && sschanlist.str().size() ) sschanlist << " ";
277  sschanlist << schan;
278  }
279  string schanlist = sschanlist.str();
280  if ( m_ShowGroups == 1 ) {
281  out << prefix << setw(wori) << iori << setw(wgrp) << igrp << ": " << schanlist << endl;
282  } else if ( m_ShowGroups == 2 ) {
283  if ( gochans.size() < igrp+1 ) gochans.resize(igrp+1, std::vector<string>(nori));
284  gochans[igrp][iori] = schanlist;
285  }
286  }
287  }
288  if ( m_ShowGroups == 2 ) {
289  for ( unsigned int igrp=0; igrp<gochans.size(); ++igrp ) {
290  for ( unsigned int iori=0; iori<nori; ++iori ) {
291  string schanlist = gochans[igrp][iori];
292  if ( schanlist.size() ) out << prefix << setw(wgrp+2) << igrp << "-" << iori << ": " << schanlist << endl;
293  }
294  }
295  }
296 
297  }
298  return out;
299 }
300 
301 //**********************************************************************
302 
304 
305 //**********************************************************************
intermediate_table::iterator iterator
int update(AdcChannelDataMap &datamap) const
#define DEFINE_ART_SERVICE_INTERFACE_IMPL(svc, iface)
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
std::string string
Definition: nybbler.cc:12
const AdcFlag AdcGood
Definition: AdcTypes.h:32
struct vector vector
Dune35tNoiseRemovalService(fhicl::ParameterSet const &pset, art::ActivityRegistry &)
intermediate_table::const_iterator const_iterator
const AdcFlag AdcStuckOn
Definition: AdcTypes.h:37
std::ostream & print(std::ostream &out=std::cout, std::string prefix="") const
float AdcSignal
Definition: AdcTypes.h:21
ROPID_t ROP
Index of the readout plane within its TPC set.
const lbne::ChannelMapService * m_pChannelMap
short AdcFlag
Definition: AdcTypes.h:29
T get(std::string const &key) const
Definition: ParameterSet.h:231
size_t size
Definition: lodepng.cpp:55
bool get_if_present(std::string const &key, T &value) const
Definition: ParameterSet.h:208
View_t View(geo::PlaneID const &pid) const
Returns the view (wire orientation) on the channels of specified TPC plane.
Class identifying a set of planes sharing readout channels.
Q_EXPORT QTSManip setw(int w)
Definition: qtextstream.h:331
AdcFilterVector signal
std::vector< std::vector< AdcChannelVector > > m_GroupChannels
std::map< AdcChannel, AdcChannelData > AdcChannelDataMap
Planes which measure W (third view for Bo, MicroBooNE, etc).
Definition: geo_types.h:127
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:402
std::vector< AdcChannel > AdcChannelVector
Definition: AdcTypes.h:51
const AdcFlag AdcStuckOff
Definition: AdcTypes.h:36
AdcSignalVector samples
unsigned int AdcChannel
Definition: AdcTypes.h:50
QTextStream & endl(QTextStream &s)
std::map< unsigned int, std::string > m_sRopChannelMap
AdcFlagVector flags