Functions
_test_noise_filtering.cxx File Reference
#include "WireCellIface/SimpleFrame.h"
#include "WireCellIface/SimpleTrace.h"
#include "WireCellSigProc/OmnibusNoiseFilter.h"
#include "WireCellSigProc/Microboone.h"
#include "WireCellSigProc/SimpleChannelNoiseDB.h"
#include "WireCellUtil/Testing.h"
#include "WireCellUtil/ExecMon.h"
#include <iostream>
#include <string>
#include <numeric>
#include "TCanvas.h"
#include "TProfile.h"
#include "TH2F.h"
#include "TFile.h"
#include "TTree.h"
#include "example-chirping-48.h"
#include "example-noisy-48.h"
#include "example-misconfig-48.h"
#include "example-rcrc-48.h"
#include "example-noisy-48-filtered.h"
#include "example-chirping-48-filtered.h"
#include "example-misconfig-48-filtered.h"
#include "example-rcrc-48-filtered.h"

Go to the source code of this file.

Functions

int main (int argc, char *argv[])
 

Function Documentation

int main ( int  argc,
char *  argv[] 
)

Definition at line 28 of file _test_noise_filtering.cxx.

29 {
30  std::vector<std::vector<float>> horigs;
31 #include "example-chirping-48.h" // run 3493 ch 720+48
32  assert(horigs.size()==48);
33  assert(horigs.at(0).size()==9594);
34 #include "example-noisy-48.h" //run 3493 ch 624+48
35  assert(horigs.size()==96);
36 #include "example-misconfig-48.h" // run3493 ch 2016+48
37  assert(horigs.size()==144);
38 #include "example-rcrc-48.h" //run3493 ch 7728+48
39  assert(horigs.size()==192);
40 
41  // std::cout << horigs.size() << " " << horigs.at(0).size() << std::endl;
42 
43  ITrace::vector traces;
44  int chindex = 624;
45  for (int ich = 0; ich!=48; ich++){
46  ITrace::ChargeSequence charges;
47  for (int itick =0; itick!=9594;itick++){
48  auto q = horigs.at(ich+48).at(itick);
49  charges.push_back(q);
50  }
51  SimpleTrace *st = new SimpleTrace(chindex+ich,0.0,charges);
52  traces.push_back(ITrace::pointer(st));
53  }
54 
55  chindex = 720;
56  for (int ich = 0; ich!=48; ich++){
57  ITrace::ChargeSequence charges;
58  for (int itick =0; itick!=9594;itick++){
59  auto q = horigs.at(ich).at(itick);
60  charges.push_back(q);
61  }
62  SimpleTrace *st = new SimpleTrace(chindex+ich,0.0,charges);
63  traces.push_back(ITrace::pointer(st));
64  }
65  chindex = 2016;
66  for (int ich = 0; ich!=48; ich++){
67  ITrace::ChargeSequence charges;
68  for (int itick =0; itick!=9594;itick++){
69  auto q = horigs.at(ich+96).at(itick);
70  charges.push_back(q);
71  }
72  SimpleTrace *st = new SimpleTrace(chindex+ich,0.0,charges);
73  traces.push_back(ITrace::pointer(st));
74  }
75  chindex = 7728;
76  for (int ich = 0; ich!=48; ich++){
77  ITrace::ChargeSequence charges;
78  for (int itick =0; itick!=9594;itick++){
79  auto q = horigs.at(ich+144).at(itick);
80  charges.push_back(q);
81  }
82  SimpleTrace *st = new SimpleTrace(chindex+ich,0.0,charges);
83  traces.push_back(ITrace::pointer(st));
84  }
85 
86 
87  SimpleFrame* sf = new SimpleFrame(0, 0, traces);
88 
89 
90 
91 
92  // S&C microboone sampling parameter database
93  const double tick = 0.5*units::microsecond;
94  const int nsamples = 9594;
95 
96  // Q&D microboone channel map
97  vector<int> uchans(2400), vchans(2400), wchans(3456);
98  const int nchans = uchans.size() + vchans.size() + wchans.size();
99  std::iota(uchans.begin(), uchans.end(), 0);
100  std::iota(vchans.begin(), vchans.end(), vchans.size());
101  std::iota(wchans.begin(), wchans.end(), vchans.size() + uchans.size());
102 
103  // Q&D nominal baseline
104  const double unombl=2048.0, vnombl=2048.0, wnombl=400.0;
105 
106  // Q&D miss-configured channel database
107  vector<int> miscfgchan;
108  const double from_gain_mVfC=4.7, to_gain_mVfC=14.0,
109  from_shaping=1.0*units::microsecond, to_shaping=2.0*units::microsecond;
110  for (int ind=2016; ind<= 2095; ++ind) { miscfgchan.push_back(ind); }
111  for (int ind=2192; ind<= 2303; ++ind) { miscfgchan.push_back(ind); }
112  for (int ind=2352; ind< 2400; ++ind) { miscfgchan.push_back(ind); }
113 
114  // hard-coded bad channels
115  vector<int> bad_channels;
116  for (int i=0;i!=wchans.size();i++){
117  if (i>=7136 - 4800 && i <=7263 - 4800){
118  if (i != 7200- 4800 && i!=7215 - 4800)
119  bad_channels.push_back(i+4800);
120  }
121  }
122 
123  // Q&D RC+RC time constant - all have same.
124  const double rcrc = 1.0*units::millisecond;
125  vector<int> rcrcchans(nchans);
126  std::iota(rcrcchans.begin(), rcrcchans.end(), 0);
127 
128  //harmonic noises
129  vector<int> harmonicchans(uchans.size() + vchans.size());
130  std::iota(harmonicchans.begin(), harmonicchans.end(), 0);
131 
132  vector<int> special_chans;
133  special_chans.push_back(2240);
134 
135  SigProc::SimpleChannelNoiseDB::mask_t h36kHz(0,169,173);
136  SigProc::SimpleChannelNoiseDB::mask_t h108kHz(0,513,516);
137  SigProc::SimpleChannelNoiseDB::mask_t hspkHz(0,17,19);
138  SigProc::SimpleChannelNoiseDB::multimask_t hharmonic;
139  hharmonic.push_back(h36kHz);
140  hharmonic.push_back(h108kHz);
141  SigProc::SimpleChannelNoiseDB::multimask_t hspecial;
142  hspecial.push_back(h36kHz);
143  hspecial.push_back(h108kHz);
144  hspecial.push_back(hspkHz);
145 
146 
147  float u_resp_array[120]={0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.364382, 0.387949, 0.411053, 0.433979, 0.456863, 0.479746, 0.502641, 0.52554, 0.548441, 0.57134, 0.591765, 0.609448, 0.626848, 0.644094, 0.661364, 0.678859, 0.695231, 0.710462, 0.726147, 0.742373, 0.761332, 0.783313, 0.806325, 0.830412, 0.857676, 0.888412, 0.920705, 0.954624, 0.990242, 1.02766, 1.06121, 1.09027, 1.12037, 1.15157, 1.18392, 1.21748, 1.25229, 1.28824, 1.32509, 1.36256, 1.40051, 1.43907, 1.47857, 1.51933, 1.56134, 1.60404, 1.72665, 1.94005, 2.16994, 2.42041, 2.69475, 3.07222, 3.67375, 4.60766, 5.91864, 7.30178, 8.3715, 8.94736, 8.93705, 8.40339, 7.2212, 5.76382, 3.8931, 1.07893, -3.52481, -11.4593, -20.4011, -29.1259, -34.9544, -36.9358, -35.3303, -31.2068, -25.8614, -20.3613, -15.3794, -11.2266, -7.96091, -5.50138, -3.71143, -2.44637, -1.57662, -0.99733, -0.62554, -0.393562, -0.249715, -0.15914, -0.100771, -0.062443, -0.037283, -0.0211508, -0.0112448, -0.00552085, -0.00245133, -0.000957821, -0.000316912, -8.51679e-05, -2.21299e-05, -1.37496e-05, -1.49806e-05, -1.36935e-05, -9.66758e-06, -5.20773e-06, -7.4787e-07, 3.71199e-06, 8.17184e-06, 1.26317e-05, 1.70916e-05, 2.15514e-05, 2.60113e-05, 3.04711e-05};
148  float v_resp_array[120]={0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0865303, 0.0925559, 0.0983619, 0.104068, 0.109739, 0.115403, 0.121068, 0.126735, 0.132403, 0.138072, 0.143739, 0.149408, 0.155085, 0.160791, 0.166565, 0.172454, 0.178514, 0.184795, 0.191341, 0.198192, 0.205382, 0.212944, 0.220905, 0.229292, 0.238129, 0.247441, 0.257256, 0.267601, 0.278502, 0.28999, 0.298745, 0.304378, 0.310105, 0.315921, 0.321818, 0.327796, 0.333852, 0.339967, 0.346098, 0.352169, 0.358103, 0.363859, 0.36945, 0.374915, 0.380261, 0.385401, 0.39016, 0.394378, 0.39804, 0.401394, 0.405145, 0.410714, 0.4205, 0.437951, 0.467841, 0.516042, 0.587738, 0.694157, 0.840763, 1.01966, 1.22894, 1.5612, 2.12348, 3.31455, 5.59355, 9.10709, 14.1756, 18.4603, 19.9517, 17.4166, 10.6683, 1.40656, -10.0638, -19.034, -23.654, -24.0558, -21.4418, -17.3229, -12.9485, -9.08912, -6.05941, -3.86946, -2.38669, -1.43678, -0.853335, -0.503951, -0.296551, -0.173029, -0.0990099, -0.0547172, -0.0287882, -0.0142758, -0.00661815, -0.00284757, -0.00115702, -0.000456456, -0.000183439, -8.04214e-05, -4.20533e-05, -2.62903e-05, -1.64098e-05, -6.68039e-06, 3.04903e-06, 1.27784e-05, 2.25079e-05, 3.22373e-05, 4.19667e-05, 5.16961e-05, 6.14255e-05, 7.11549e-05};
149  WireCell::Waveform::realseq_t u_resp(nsamples);
150  WireCell::Waveform::realseq_t v_resp(nsamples);
151  for (int i=0;i!=120;i++){
152  u_resp.at(i) = u_resp_array[i];
153  v_resp.at(i) = v_resp_array[i];
154  }
155  WireCell::Waveform::compseq_t u_resp_freq = WireCell::Waveform::dft(u_resp);
156  WireCell::Waveform::compseq_t v_resp_freq = WireCell::Waveform::dft(v_resp);
157 
158 
159  int uplane_time_shift = 79;
160  int vplane_time_shift = 82;
161 
162 
163 
164  // do the coherent subtraction
165 
166  std::vector< std::vector<int> > channel_groups;
167  for (int i=0;i!=172;i++){
168  //for (int i=150;i!=151;i++){
169  std::vector<int> channel_group;
170  for (int j=0;j!=48;j++){
171  channel_group.push_back(i*48+j);
172  }
173  channel_groups.push_back(channel_group);
174  }
175 
176  // Load up components. Note, in a real app this is done as part
177  // of factory + configurable and driven by user configuration.
178 
179  auto noise = new SigProc::SimpleChannelNoiseDB;
180  // initialize
181  noise->set_sampling(tick, nsamples);
182  // set nominal baseline
183  noise->set_nominal_baseline(uchans, unombl);
184  noise->set_nominal_baseline(vchans, vnombl);
185  noise->set_nominal_baseline(wchans, wnombl);
186 
187  noise->set_response(uchans,u_resp_freq);
188  noise->set_response(vchans,v_resp_freq);
189 
190  noise->set_response_offset(uchans,uplane_time_shift);
191  noise->set_response_offset(vchans,vplane_time_shift);
192 
193  noise->set_pad_window_front(uchans,20);
194  //noise->set_pad_window_front(uchans,10);
195  noise->set_pad_window_back(uchans,10);
196  noise->set_pad_window_front(vchans,10);
197  noise->set_pad_window_back(vchans,10);
198  noise->set_pad_window_front(wchans,10);
199  noise->set_pad_window_back(wchans,10);
200 
201 
202  // set misconfigured channels
203  noise->set_gains_shapings(miscfgchan, from_gain_mVfC, to_gain_mVfC, from_shaping, to_shaping);
204  // do the RCRC
205  noise->set_rcrc_constant(rcrcchans, rcrc);
206 
207 
208  // fill in the rms cut ...
209  // before Hardware Fix
210  for (int i=0;i!=uchans.size();i++){
211  if (uchans.at(i)<100){
212  noise->set_min_rms_cut_one(uchans.at(i),1);
213  noise->set_max_rms_cut_one(uchans.at(i),5);
214  }else if (uchans.at(i)<2000){
215  noise->set_min_rms_cut_one(uchans.at(i),1.9);
216  noise->set_max_rms_cut_one(uchans.at(i),11);
217  }else{
218  noise->set_min_rms_cut_one(uchans.at(i),0.9);
219  noise->set_max_rms_cut_one(uchans.at(i),5);
220  }
221  }
222  for (int i=0;i!=vchans.size();i++){
223  if (vchans.at(i)<290+uchans.size()){
224  noise->set_min_rms_cut_one(vchans.at(i),1);
225  noise->set_max_rms_cut_one(vchans.at(i),5);
226  }else if (vchans.at(i)<2200+uchans.size()){
227  noise->set_min_rms_cut_one(vchans.at(i),1.9);
228  noise->set_max_rms_cut_one(vchans.at(i),11);
229  }else{
230  noise->set_min_rms_cut_one(vchans.at(i),1);
231  noise->set_max_rms_cut_one(vchans.at(i),5);
232  }
233  }
234  noise->set_min_rms_cut(wchans,1.3);
235  noise->set_max_rms_cut(wchans,8.0);
236 
237 
238 
239 
240 
241 
242 
243  // set initial bad channels
244  noise->set_bad_channels(bad_channels);
245  // set the harmonic filter
246  noise->set_filter(harmonicchans,hharmonic);
247  noise->set_filter(special_chans,hspecial);
248  noise->set_channel_groups(channel_groups);
249 
250  shared_ptr<WireCell::IChannelNoiseDatabase> noise_sp(noise);
251 
252  auto one = new SigProc::Microboone::OneChannelNoise;
253  one->set_channel_noisedb(noise_sp);
254  shared_ptr<WireCell::IChannelFilter> one_sp(one);
255 
256  auto many = new SigProc::Microboone::CoherentNoiseSub;
257  many->set_channel_noisedb(noise_sp);
258  shared_ptr<WireCell::IChannelFilter> many_sp(many);
259 
260 
261  SigProc::OmnibusNoiseFilter bus;
262  bus.set_channel_filters({one_sp});
263  bus.set_grouped_filters({many_sp});
264  bus.set_channel_noisedb(noise_sp);
265 
266  IFrame::pointer frame = IFrame::pointer(sf);
267 
268  IFrame::pointer quiet;
269  bus(frame, quiet);
270  Assert(quiet);
271 
272  std::vector<std::vector<float>> hfilts;
273 
274 #include "example-noisy-48-filtered.h"
275  assert(hfilts.size()==48);
276  assert(hfilts.at(0).size()==9594);
277 #include "example-chirping-48-filtered.h"
278  assert(hfilts.size()==96);
279 #include "example-misconfig-48-filtered.h"
280  assert(hfilts.size()==144);
281 #include "example-rcrc-48-filtered.h"
282  assert(hfilts.size()==192);
283  // test ...
284  auto traces1 = quiet->traces();
285  int counter = 0;
286 
287  for (auto trace : *traces1.get()) {
288  int tbin = trace->tbin();
289  int ch = trace->channel();
290  auto charges = trace->charge();
291  // std::cout << ch << " " << counter << " " << charges.size() << " " << hfilts.at(counter).size() << " " << charges.at(0) << " " << hfilts.at(counter).at(0) << std::endl;
292  for (int i=0;i!=9594;i++){
293  const double diff = fabs(charges.at(i) - hfilts.at(counter).at(i));
294  if (diff >= 1.0) {
295  std::cerr << "Warning: normally would assert: "
296  << "counter: " << counter << " ch:" << i
297  << " diff: " << diff
298  << " charge:" << charges.at(i)
299  << " hfilts:" << hfilts.at(counter).at(i)
300  << std::endl;
301  }
302  //Assert( diff < 1.0 );
303  }
304 
305 
306  counter ++;
307  }
308 
309  return 0;
310 }
WireCell::Waveform::realseq_t
Sequence< real_t > realseq_t
Definition: Waveform.h:31
WireCell::SigProc::OmnibusNoiseFilter::set_channel_filters
void set_channel_filters(std::vector< WireCell::IChannelFilter::pointer > filters)
Definition: OmnibusNoiseFilter.h:52
util::quantities::microsecond
microsecond_as<> microsecond
Type of time stored in microseconds, in double precision.
Definition: spacetime.h:114
vector
struct vector vector
WireCell::SigProc::Microboone::ConfigFilterBase::set_channel_noisedb
void set_channel_noisedb(WireCell::IChannelNoiseDatabase::pointer ndb)
Definition: Microboone.h:49
tick
const double tick
Definition: test_binneddiffusion.cxx:51
util::quantities::millisecond
millisecond_as<> millisecond
Type of time stored in milliseconds, in double precision.
Definition: spacetime.h:97
Assert
#define Assert
Definition: Testing.h:7
WireCell::Waveform::dft
compseq_t dft(realseq_t seq)
Definition: Waveform.cxx:141
util::counter
auto counter(T begin, T end)
Returns an object to iterate values from begin to end in a range-for loop.
Definition: counter.h:285
WireCell::SigProc::OmnibusNoiseFilter::set_channel_noisedb
void set_channel_noisedb(WireCell::IChannelNoiseDatabase::pointer ndb)
Definition: OmnibusNoiseFilter.h:61
WireCell::SigProc::OmnibusNoiseFilter::set_grouped_filters
void set_grouped_filters(std::vector< WireCell::IChannelFilter::pointer > filters)
Definition: OmnibusNoiseFilter.h:58
pointer
const GenericPointer< typename T::ValueType > & pointer
Definition: pointer.h:1124
WireCell::Waveform::compseq_t
Sequence< complex_t > compseq_t
A complex-valued sequence, eg for discrete spectrum powers.
Definition: Waveform.h:34
WireCell::ITrace::ChargeSequence
std::vector< float > ChargeSequence
Sequential collection of charge.
Definition: ITrace.h:21
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
WireCell::SigProc::SimpleChannelNoiseDB::set_sampling
void set_sampling(double tick=0.5 *units::us, int nsamples=9600)
Definition: SimpleChannelNoiseDB.cxx:213
WireCell::SigProc::SimpleChannelNoiseDB::mask_t
std::tuple< double, int, int > mask_t
Definition: SimpleChannelNoiseDB.h:119
Generated by   doxygen 1.8.11