Classes | Functions
WireCell::SigProc::Microboone Namespace Reference

Classes

class  ADCBitShift
 
class  CoherentNoiseSub
 
class  ConfigFilterBase
 
class  OneChannelNoise
 
class  OneChannelStatus
 

Functions

bool Chirp_raise_baseline (WireCell::Waveform::realseq_t &sig, int bin1, int bin2)
 
bool SignalFilter (WireCell::Waveform::realseq_t &sig)
 
float CalcRMSWithFlags (const WireCell::Waveform::realseq_t &sig)
 
bool RawAdapativeBaselineAlg (WireCell::Waveform::realseq_t &sig)
 
bool RemoveFilterFlags (WireCell::Waveform::realseq_t &sig)
 
bool NoisyFilterAlg (WireCell::Waveform::realseq_t &spec, float min_rms, float max_rms)
 
std::vector< std::vector< int > > SignalProtection (WireCell::Waveform::realseq_t &sig, const WireCell::Waveform::compseq_t &respec, int res_offset, int pad_f, int pad_b, float upper_decon_limit=0.02, float decon_lf_cutoff=0.08, float upper_adc_limit=15, float protection_factor=5.0, float min_adc_limit=50)
 
bool Subtract_WScaling (WireCell::IChannelFilter::channel_signals_t &chansig, const WireCell::Waveform::realseq_t &medians, const WireCell::Waveform::compseq_t &respec, int res_offset, std::vector< std::vector< int > > &rois, float upper_decon_limit1=0.08, float roi_min_max_ratio=0.8)
 

Function Documentation

float WireCell::SigProc::Microboone::CalcRMSWithFlags ( const WireCell::Waveform::realseq_t &  sig)

Definition at line 625 of file Microboone.cxx.

626 {
627  float theRMS = 0.0;
628  //int waveformSize = sig.size();
629 
630  WireCell::Waveform::realseq_t temp;
631  for (size_t i=0;i!=sig.size();i++){
632  if (sig.at(i) < 4096) temp.push_back(sig.at(i));
633  }
634  float par[3];
635  if (temp.size()>0) {
636  par[0] = WireCell::Waveform::percentile_binned(temp,0.5 - 0.34);
637  par[1] = WireCell::Waveform::percentile_binned(temp,0.5);
638  par[2] = WireCell::Waveform::percentile_binned(temp,0.5 + 0.34);
639 
640  // std::cout << par[0] << " " << par[1] << " " << par[2] << std::endl;
641 
642  theRMS = sqrt((pow(par[2]-par[1],2)+pow(par[1]-par[0],2))/2.);
643  }
644 
645  return theRMS;
646 }
WireCell::Waveform::realseq_t
Sequence< real_t > realseq_t
Definition: Waveform.h:31
cet::pow
constexpr T pow(T x)
Definition: pow.h:72
WireCell::Waveform::percentile_binned
real_t percentile_binned(realseq_t &wave, real_t percentage)
Definition: Waveform.cxx:93
bool WireCell::SigProc::Microboone::Chirp_raise_baseline ( WireCell::Waveform::realseq_t &  sig,
int  bin1,
int  bin2 
)

Definition at line 615 of file Microboone.cxx.

616 {
617  if (bin1 < 0 ) bin1 = 0;
618  if (bin2 > (int)sig.size()) bin2 = sig.size();
619  for (int i=bin1; i<bin2;i++) {
620  sig.at(i) = 10000.0;
621  }
622  return true;
623 }
bool WireCell::SigProc::Microboone::NoisyFilterAlg ( WireCell::Waveform::realseq_t &  spec,
float  min_rms,
float  max_rms 
)

Definition at line 542 of file Microboone.cxx.

543 {
544  const double rmsVal = Microboone::CalcRMSWithFlags(sig);
545 
546  //std::cout << rmsVal << std::endl;
547 
548  // int planeNum,channel_no;
549  // if (ch < 2400) {
550  // planeNum = 0;
551  // channel_no = ch;
552  // }
553  // else if (ch < 4800) {
554  // planeNum = 1;
555  // channel_no = ch - 2400;
556  // }
557  // else {
558  // planeNum = 2;
559  // channel_no = ch - 4800;
560  // }
561 
562  // //std::cout << rmsVal << " " << ch << std::endl;
563 
564  // double maxRMSCut[3] = {10.0,10.0,5.0};
565  // double minRMSCut[3] = {2,2,2};
566 
567  // if (planeNum == 0) {
568  // if (channel_no < 100) {
569  // maxRMSCut[0] = 5;
570  // minRMSCut[0] = 1;
571  // }
572  // else if (channel_no >= 100 && channel_no<2000) {
573  // maxRMSCut[0] = 11; // increase the threshold slightly ...
574  // minRMSCut[0] = 1.9;
575  // }
576  // else if (channel_no >= 2000 && channel_no < 2400) {
577  // maxRMSCut[0] = 5;
578  // minRMSCut[0] = 0.9; // take into account FT-1 channel ...
579  // }
580  // }
581  // else if (planeNum == 1) {
582  // if (channel_no <290){
583  // maxRMSCut[1] = 5;
584  // minRMSCut[1] = 1;
585  // }
586  // else if (channel_no>=290 && channel_no < 2200) {
587  // maxRMSCut[1] = 11;
588  // minRMSCut[1] = 1.9;
589  // }
590  // else if (channel_no >=2200) {
591  // maxRMSCut[1] = 5;
592  // minRMSCut[1] = 1;
593  // }
594  // }
595  // else if (planeNum == 2) {
596  // maxRMSCut[2] = 8;
597  // minRMSCut[2] = 1.3; // reduce threshold to take into account the adaptive baseline ...
598  // }
599 
600  //if(rmsVal > maxRMSCut[planeNum] || rmsVal < minRMSCut[planeNum]) {
601  if(rmsVal > max_rms || rmsVal < min_rms) {
602  int numBins = sig.size();
603  for(int i = 0; i < numBins; i++) {
604  sig.at(i) = 10000.0;
605  }
606 
607  return true;
608  }
609 
610  return false;
611 }
WireCell::SigProc::Microboone::CalcRMSWithFlags
float CalcRMSWithFlags(const WireCell::Waveform::realseq_t &sig)
Definition: Microboone.cxx:625
bool WireCell::SigProc::Microboone::RawAdapativeBaselineAlg ( WireCell::Waveform::realseq_t &  sig)

Definition at line 693 of file Microboone.cxx.

694 {
695  const int windowSize = 20;
696  const int numBins = sig.size();
697  int minWindowBins = windowSize/2;
698 
699  std::vector<double> baselineVec(numBins, 0.0);
700  std::vector<bool> isFilledVec(numBins, false);
701 
702  int numFlaggedBins = 0;
703 
704  for(int j = 0; j < numBins; j++) {
705  if(sig.at(j) == 10000.0) {
706  numFlaggedBins++;
707  }
708  }
709  if(numFlaggedBins == numBins) {
710  return true; // Eventually replace this with flag check
711  }
712 
713  double baselineVal = 0.0;
714  int windowBins = 0;
715  //int index;
716  double ADCval = 0.0;
717  for(int j = 0; j <= windowSize/2; j++) {
718  ADCval = sig.at(j);
719  if(ADCval < 4096.0) {
720  baselineVal += ADCval;
721  windowBins++;
722  }
723  }
724 
725  if(windowBins == 0) {
726  baselineVec[0] = 0.0;
727  }
728  else {
729  baselineVec[0] = baselineVal/((double) windowBins);
730  }
731 
732  if(windowBins < minWindowBins) {
733  isFilledVec[0] = false;
734  }
735  else {
736  isFilledVec[0] = true;
737  }
738 
739  for(int j = 1; j < numBins; j++) {
740  int oldIndex = j-windowSize/2-1;
741  int newIndex = j+windowSize/2;
742 
743  if(oldIndex >= 0) {
744  ADCval = sig.at(oldIndex);
745  if(ADCval < 4096.0) {
746  baselineVal -= sig.at(oldIndex);
747  windowBins--;
748  }
749  }
750  if(newIndex < numBins) {
751  ADCval = sig.at(newIndex);
752  if(ADCval < 4096) {
753  baselineVal += sig.at(newIndex);
754  windowBins++;
755  }
756  }
757 
758  if(windowBins == 0) {
759  baselineVec[j] = 0.0;
760  }
761  else {
762  baselineVec[j] = baselineVal/windowBins;
763  }
764 
765  if(windowBins < minWindowBins) {
766  isFilledVec[j] = false;
767  }
768  else {
769  isFilledVec[j] = true;
770  }
771  }
772 
773  for(int j = 0; j < numBins; j++) {
774  bool downFlag = false;
775  bool upFlag = false;
776 
777  ADCval = sig.at(j);
778  if(ADCval != 10000.0) {
779  if(isFilledVec[j] == false) {
780  int downIndex = j;
781  while((isFilledVec[downIndex] == false) && (downIndex > 0) && (sig.at(downIndex) != 10000.0)) {
782  downIndex--;
783  }
784 
785  if(isFilledVec[downIndex] == false) {
786  downFlag = true;
787  }
788 
789  int upIndex = j;
790  while((isFilledVec[upIndex] == false) && (upIndex < numBins-1) && (sig.at(upIndex) != 10000.0)) {
791  upIndex++;
792  }
793 
794  if(isFilledVec[upIndex] == false) {
795  upFlag = true;
796  }
797 
798  if((downFlag == false) && (upFlag == false)) {
799  baselineVec[j] = ((j-downIndex)*baselineVec[downIndex]+(upIndex-j)*baselineVec[upIndex])/((double) upIndex-downIndex);
800  }
801  else if((downFlag == true) && (upFlag == false)) {
802  baselineVec[j] = baselineVec[upIndex];
803  }
804  else if((downFlag == false) && (upFlag == true)) {
805  baselineVec[j] = baselineVec[downIndex];
806  }
807  else {
808  baselineVec[j] = 0.0;
809  }
810  }
811 
812  sig.at(j) = ADCval -baselineVec[j];
813  }
814  }
815 
816  return true;
817 
818 }
bool WireCell::SigProc::Microboone::RemoveFilterFlags ( WireCell::Waveform::realseq_t &  sig)

Definition at line 821 of file Microboone.cxx.

822 {
823  int numBins = sig.size();
824  for(int i = 0; i < numBins; i++) {
825  double ADCval = sig.at(i);
826  if (ADCval > 4096.0) {
827  if(ADCval > 10000.0) {
828  sig.at(i) = ADCval-20000.0;
829  }
830  else {
831  sig.at(i) = 0.0;
832  }
833 
834  }
835  }
836 
837  return true;
838 }
bool WireCell::SigProc::Microboone::SignalFilter ( WireCell::Waveform::realseq_t &  sig)

Definition at line 648 of file Microboone.cxx.

649 {
650  const double sigFactor = 4.0;
651  const int padBins = 8;
652 
653  float rmsVal = Microboone::CalcRMSWithFlags(sig);
654  float sigThreshold = sigFactor*rmsVal;
655 
656  float ADCval;
657  std::vector<bool> signalRegions;
658  int numBins = sig.size();
659 
660  for (int i = 0; i < numBins; i++) {
661  ADCval = sig.at(i);
662  if (((ADCval > sigThreshold) || (ADCval < -1.0*sigThreshold)) && (ADCval < 4096.0)) {
663  signalRegions.push_back(true);
664  }
665  else {
666  signalRegions.push_back(false);
667  }
668  }
669 
670  for(int i = 0; i < numBins; i++) {
671  if(signalRegions[i] == true) {
672  int bin1 = i - padBins;
673  if (bin1 < 0 ) {
674  bin1 = 0;
675  }
676  int bin2 = i + padBins;
677  if (bin2 > numBins) {
678  bin2 = numBins;
679  }
680 
681  for(int j = bin1; j < bin2; j++) {
682  ADCval = sig.at(j);
683  if(ADCval < 4096.0) {
684  sig.at(j) = sig.at(j)+20000.0;
685  }
686  }
687  }
688  }
689  return true;
690 }
WireCell::SigProc::Microboone::CalcRMSWithFlags
float CalcRMSWithFlags(const WireCell::Waveform::realseq_t &sig)
Definition: Microboone.cxx:625
std::vector< std::vector< int > > WireCell::SigProc::Microboone::SignalProtection ( WireCell::Waveform::realseq_t &  sig,
const WireCell::Waveform::compseq_t &  respec,
int  res_offset,
int  pad_f,
int  pad_b,
float  upper_decon_limit = 0.02,
float  decon_lf_cutoff = 0.08,
float  upper_adc_limit = 15,
float  protection_factor = 5.0,
float  min_adc_limit = 50 
)

Definition at line 232 of file Microboone.cxx.

233 {
234 
235 
236  // WireCell::Waveform::realseq_t temp1;
237  // for (int i=0;i!=medians.size();i++){
238  // if (fabs(medians.at(i) - mean) < 4.5*rms)
239  // temp1.push_back(medians.at(i));
240  // }
241  // temp = WireCell::Waveform::mean_rms(temp1);
242  // mean = temp.first;
243  // rms = temp.second;
244 
245  //std::cout << temp.first << " " << temp.second << std::endl;
246  const int nbin = medians.size();
247 
248  // const int protection_factor = 5.0;
249  // move to input ...
250  // const float upper_decon_limit = 0.05;
251  // const float upper_adc_limit = 15;
252  //const float min_adc_limit = 50;
253 
254 
255  std::vector<bool> signalsBool;
256  signalsBool.resize(nbin, false);
257 
258 
259 
260 
261  // calculate the RMS
262  std::pair<double,double> temp = Derivations::CalcRMS(medians);
263  double mean = temp.first;
264  double rms = temp.second;
265 
266  // std::cout << mean << " " << rms << " " << respec.size() << " " << res_offset << " " << pad_f << " " << pad_b << " " << respec.at(0) << std::endl;
267 
268  float limit;
269  if (protection_factor*rms > upper_adc_limit){
270  limit = protection_factor*rms;
271  }else{
272  limit = upper_adc_limit;
273  }
274  if (min_adc_limit < limit){
275  limit = min_adc_limit;
276  }
277 
278  // std::cout << "Xin " << protection_factor << " " << mean << " " << rms *protection_factor << " " << upper_adc_limit << " " << decon_lf_cutoff << " " << upper_decon_limit << std::endl;
279 
280  for (int j=0;j!=nbin;j++) {
281  float content = medians.at(j);
282  if (fabs(content-mean)>limit){
283  //protection_factor*rms) {
284  // medians.at(j) = 0;
285  signalsBool.at(j) = true;
286  // add the front and back padding
287  for (int k=0;k!=pad_b;k++){
288  int bin = j+k+1;
289  if (bin > nbin-1) bin = nbin-1;
290  signalsBool.at(bin) = true;
291  }
292  for (int k=0;k!=pad_f;k++){
293  int bin = j-k-1;
294  if (bin <0) { bin = 0; }
295  signalsBool.at(bin) = true;
296  }
297  }
298  }
299 
300 
301  //std::cout << "Xin: " << respec.size() << " " << res_offset << std::endl;
302  // the deconvolution protection code ...
303  if (respec.size() > 0 && (respec.at(0).real()!=1 || respec.at(0).imag()!=0) && res_offset!=0){
304  //std::cout << nbin << std::endl;
305 
306  WireCell::Waveform::compseq_t medians_freq = WireCell::Waveform::dft(medians);
307  WireCell::Waveform::shrink(medians_freq,respec);
308 
309  for (size_t i=0;i!=medians_freq.size();i++){
310  double freq;
311  // assuming 2 MHz digitization
312  if (i <medians_freq.size()/2.){
313  freq = i/(1.*medians_freq.size())*2.;
314  }else{
315  freq = (medians_freq.size() - i)/(1.*medians_freq.size())*2.;
316  }
317  std::complex<float> factor = filter_time(freq)*filter_low(freq, decon_lf_cutoff);
318  medians_freq.at(i) = medians_freq.at(i) * factor;
319  }
320  WireCell::Waveform::realseq_t medians_decon = WireCell::Waveform::idft(medians_freq);
321 
322  temp = Derivations::CalcRMS(medians_decon);
323  mean = temp.first;
324  rms = temp.second;
325 
326  if (protection_factor*rms > upper_decon_limit){
327  limit = protection_factor*rms;
328  }else{
329  limit = upper_decon_limit;
330  }
331 
332 
333  // std::cout << "Xin: " << protection_factor << " " << rms << " " << upper_decon_limit << std::endl;
334 
335 
336 
337  for (int j=0;j!=nbin;j++) {
338  float content = medians_decon.at(j);
339  if ((content-mean)>limit){
340  int time_bin = j + res_offset;
341  if (time_bin >= nbin) time_bin -= nbin;
342  // medians.at(time_bin) = 0;
343  signalsBool.at(time_bin) = true;
344  // add the front and back padding
345  for (int k=0;k!=pad_b;k++){
346  int bin = time_bin+k+1;
347  if (bin > nbin-1) bin = nbin-1;
348  signalsBool.at(bin) = true;
349  }
350  for (int k=0;k!=pad_f;k++){
351  int bin = time_bin-k-1;
352  if (bin <0) { bin = 0; }
353  signalsBool.at(bin) = true;
354  }
355  }
356  }
357 
358 
359  // // second-level decon ...
360  // medians_freq = WireCell::Waveform::dft(medians);
361  // WireCell::Waveform::realseq_t respec_time = WireCell::Waveform::idft(respec);
362  // for (size_t i=0;i!=respec_time.size();i++){
363  // if (respec_time.at(i)<0) respec_time.at(i) = 0;
364  // }
365  // WireCell::Waveform::compseq_t respec_freq = WireCell::Waveform::dft(respec_time);
366  // WireCell::Waveform::shrink(medians_freq,respec_freq);
367  // for (size_t i=0;i!=medians_freq.size();i++){
368  // double freq;
369  // // assuming 2 MHz digitization
370  // if (i <medians_freq.size()/2.){
371  // freq = i/(1.*medians_freq.size())*2.;
372  // }else{
373  // freq = (medians_freq.size() - i)/(1.*medians_freq.size())*2.;
374  // }
375  // std::complex<float> factor = filter_time(freq)*filter_low(freq, decon_lf_cutoff);
376  // medians_freq.at(i) = medians_freq.at(i) * factor;
377  // }
378  // medians_decon = WireCell::Waveform::idft(medians_freq);
379 
380  // temp = Derivations::CalcRMS(medians_decon);
381  // mean = temp.first;
382  // rms = temp.second;
383 
384  // // if (protection_factor*rms > upper_decon_limit){
385  // limit = protection_factor*rms;
386  // // }else{
387  // // limit = upper_decon_limit;
388  // // }
389 
390  // for (int j=0;j!=nbin;j++) {
391  // float content = medians_decon.at(j);
392  // if ((content-mean)>limit){
393  // int time_bin = j + res_offset;
394  // if (time_bin >= nbin) time_bin -= nbin;
395  // // medians.at(time_bin) = 0;
396  // signalsBool.at(time_bin) = true;
397  // // add the front and back padding
398  // for (int k=0;k!=pad_b;k++){
399  // int bin = time_bin+k+1;
400  // if (bin > nbin-1) bin = nbin-1;
401  // signalsBool.at(bin) = true;
402  // }
403  // for (int k=0;k!=pad_f;k++){
404  // int bin = time_bin-k-1;
405  // if (bin <0) { bin = 0; }
406  // signalsBool.at(bin) = true;
407  // }
408  // }
409  // }
410 
411 
412  }
413 
414  // {
415  // partition waveform indices into consecutive regions with
416  // signalsBool true.
417  std::vector< std::vector<int> > rois;
418  bool inside = false;
419  for (int ind=0; ind<nbin; ++ind) {
420  if (inside) {
421  if (signalsBool[ind]) { // still inside
422  rois.back().push_back(ind);
423  }else{
424  inside = false;
425  }
426  }
427  else { // outside the Rio
428  if (signalsBool[ind]) { // just entered ROI
429  std::vector<int> roi;
430  roi.push_back(ind);
431  rois.push_back(roi);
432  inside = true;
433  }
434  }
435  }
436 
437 
438 
439  std::map<int, bool> flag_replace;
440  for (auto roi: rois){
441  flag_replace[roi.front()] = true;
442  }
443 
444  if (respec.size() > 0 && (respec.at(0).real()!=1 || respec.at(0).imag()!=0) && res_offset!=0){
445  for (auto roi: rois){
446  flag_replace[roi.front()] = false;
447  }
448  }
449 
450 
451  // // use ROI to get a new waveform
452  // WireCell::Waveform::realseq_t medians_roi(nbin,0);
453  // for (auto roi: rois){
454  // const int bin0 = std::max(roi.front()-1, 0);
455  // const int binf = std::min(roi.back()+1, nbin-1);
456  // const double m0 = medians[bin0];
457  // const double mf = medians[binf];
458  // const double roi_run = binf - bin0;
459  // const double roi_rise = mf - m0;
460  // for (auto bin : roi) {
461  // const double m = m0 + (bin - bin0)/roi_run*roi_rise;
462  // medians_roi.at(bin) = medians.at(bin) - m;
463  // }
464  // }
465  // // do the deconvolution with a very loose low-frequency filter
466  // WireCell::Waveform::compseq_t medians_roi_freq = WireCell::Waveform::dft(medians_roi);
467  // WireCell::Waveform::shrink(medians_roi_freq,respec);
468  // for (size_t i=0;i!=medians_roi_freq.size();i++){
469  // double freq;
470  // // assuming 2 MHz digitization
471  // if (i <medians_roi_freq.size()/2.){
472  // freq = i/(1.*medians_roi_freq.size())*2.;
473  // }else{
474  // freq = (medians_roi_freq.size() - i)/(1.*medians_roi_freq.size())*2.;
475  // }
476  // std::complex<float> factor = filter_time(freq)*filter_low_loose(freq);
477  // medians_roi_freq.at(i) = medians_roi_freq.at(i) * factor;
478  // }
479  // WireCell::Waveform::realseq_t medians_roi_decon = WireCell::Waveform::idft(medians_roi_freq);
480 
481  // // judge if a roi is good or not ...
482  // //shift things back properly
483  // for (auto roi: rois){
484  // const int bin0 = std::max(roi.front()-1, 0);
485  // const int binf = std::min(roi.back()+1, nbin-1);
486  // flag_replace[roi.front()] = false;
487 
488  // double max_val = 0;
489  // double min_val = 0;
490  // // double max_adc_val=0;
491  // // double min_adc_val=0;
492 
493  // for (int i=bin0; i<=binf; i++){
494  // int time_bin = i-res_offset;
495  // if (time_bin <0) time_bin += nbin;
496  // if (time_bin >=nbin) time_bin -= nbin;
497 
498  // if (i==bin0){
499  // max_val = medians_roi_decon.at(time_bin);
500  // min_val = medians_roi_decon.at(time_bin);
501  // // max_adc_val = medians.at(i);
502  // // min_adc_val = medians.at(i);
503  // }else{
504  // if (medians_roi_decon.at(time_bin) > max_val) max_val = medians_roi_decon.at(time_bin);
505  // if (medians_roi_decon.at(time_bin) < min_val) min_val = medians_roi_decon.at(time_bin);
506  // // if (medians.at(i) > max_adc_val) max_adc_val = medians.at(i);
507  // // if (medians.at(i) < min_adc_val) min_adc_val = medians.at(i);
508  // }
509  // }
510 
511  // //std::cout << "Xin: " << upper_decon_limit1 << std::endl;
512  // // if ( max_val > upper_decon_limit1)
513  // // if ( max_val > 0.04 && fabs(min_val) < 0.6*max_val)
514  // //if (max_val > 0.06 && fabs(min_val) < 0.6*max_val)
515  // if (max_val > 0.06)
516  // flag_replace[roi.front()] = true;
517  // }
518  // }
519 
520  // Replace medians for above regions with interpolation on values
521  // just outside each region.
522  for (auto roi : rois) {
523  // original code used the bins just outside the ROI
524  const int bin0 = std::max(roi.front()-1, 0);
525  const int binf = std::min(roi.back()+1, nbin-1);
526  if (flag_replace[roi.front()]){
527  const double m0 = medians[bin0];
528  const double mf = medians[binf];
529  const double roi_run = binf - bin0;
530  const double roi_rise = mf - m0;
531  for (auto bin : roi) {
532  const double m = m0 + (bin - bin0)/roi_run*roi_rise;
533  medians.at(bin) = m;
534  }
535  }
536  }
537 
538 
539  return rois;
540 }
genie::units::m
static const double m
Definition: Units.h:79
WireCell::Waveform::realseq_t
Sequence< real_t > realseq_t
Definition: Waveform.h:31
cet::sqlite::rms
double rms(sqlite3 *db, std::string const &table_name, std::string const &column_name)
Definition: statistics.cc:39
WireCell::SigProc::Derivations::CalcRMS
std::pair< double, double > CalcRMS(const WireCell::Waveform::realseq_t &signal)
Definition: Derivations.cxx:7
WireCell::Waveform::dft
compseq_t dft(realseq_t seq)
Definition: Waveform.cxx:141
WireCell::Waveform::shrink
void shrink(Sequence< Val > &seq, const Sequence< Val > &other)
Shrink (divide) seq values by values from the other sequence.
Definition: Waveform.h:162
max
static int max(int a, int b)
Definition: fortranscanner.cpp:60366
WireCell::Waveform::idft
realseq_t idft(compseq_t spec)
Definition: Waveform.cxx:149
cet::sqlite::min
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
bin
QTextStream & bin(QTextStream &s)
Definition: qtextstream.cpp:2006
filter_low
double filter_low(double freq, double cut_off=0.08)
Definition: Microboone.cxx:41
WireCell::Waveform::compseq_t
Sequence< complex_t > compseq_t
A complex-valued sequence, eg for discrete spectrum powers.
Definition: Waveform.h:34
cet::sqlite::mean
double mean(sqlite3 *db, std::string const &table_name, std::string const &column_name)
Definition: statistics.cc:15
filter_time
double filter_time(double freq)
Definition: Microboone.cxx:35
bool WireCell::SigProc::Microboone::Subtract_WScaling ( WireCell::IChannelFilter::channel_signals_t &  chansig,
const WireCell::Waveform::realseq_t &  medians,
const WireCell::Waveform::compseq_t &  respec,
int  res_offset,
std::vector< std::vector< int > > &  rois,
float  upper_decon_limit1 = 0.08,
float  roi_min_max_ratio = 0.8 
)

Definition at line 55 of file Microboone.cxx.

62 {
63  double ave_coef = 0;
64  double_t ave_coef1 = 0;
65  std::map<int,double> coef_all;
66 
67  const int nbin = medians.size();
68 
69  for (auto it: chansig) {
70  int ch = it.first;
71  WireCell::IChannelFilter::signal_t& signal = it.second;
72 
73 
74  double sum2 = 0;
75  double sum3 = 0;
76  double coef = 0;
77 
78  std::pair<double,double> temp = Derivations::CalcRMS(signal);
79  //std::pair<double,double> temp = WireCell::Waveform::mean_rms(signal);
80 
81  // if ( abs(ch-6117)<5)
82  // std::cout << ch << " " << temp.first << " " << temp.second << " " << std::endl;
83 
84  for (int j=0;j!=nbin;j++){
85  if (fabs(signal.at(j)) < 4 * temp.second){
86  sum2 += signal.at(j) * medians.at(j);
87  sum3 += medians.at(j) * medians.at(j);
88  }
89  }
90  if (sum3 >0) {
91  coef = sum2/sum3;
92  }
93  // protect against the extreme cases
94  // if (coef < 0.6) coef = 0.6;
95  // if (coef > 1.5) coef = 1.5;
96 
97  coef_all[ch] = coef;
98  if (coef != 0){ // FIXME: expect some fluctuation?
99  ave_coef += coef;
100  ave_coef1 ++;
101  }
102  }
103  if (ave_coef1 > 0) {
104  ave_coef = ave_coef / ave_coef1;
105  }
106 
107  for (auto it: chansig) {
108  int ch = it.first;
109  WireCell::IChannelFilter::signal_t& signal = it.second;
110  float scaling;
111  if (ave_coef != 0 ){
112  scaling = coef_all[ch]/ave_coef;
113  // add some protections ...
114  if (scaling < 0) scaling = 0;
115  // if (scaling < 0.5 && scaling > 0.3) scaling = 0.5;
116  if (scaling > 1.5) scaling = 1.5;
117  }else{
118  scaling = 0;
119  }
120  // if ( abs(ch-6117)<5)
121  // std::cout << ch << " " << scaling << " " << std::endl;
122  //scaling = 1.0;
123 
124 
125  if (respec.size() > 0 && (respec.at(0).real()!=1 || respec.at(0).imag()!=0) && res_offset!=0){
126  int nbin = signal.size();
127  WireCell::Waveform::realseq_t signal_roi(nbin,0);
128  for (auto roi: rois){
129  const int bin0 = std::max(roi.front()-1, 0);
130  const int binf = std::min(roi.back()+1, nbin-1);
131  const double m0 = signal[bin0];
132  const double mf = signal[binf];
133  const double roi_run = binf - bin0;
134  const double roi_rise = mf - m0;
135  for (auto bin : roi) {
136  const double m = m0 + (bin - bin0)/roi_run*roi_rise;
137  signal_roi.at(bin) = signal.at(bin) - m;
138  }
139  }
140 
141  // do the deconvolution with a very loose low-frequency filter
142  WireCell::Waveform::compseq_t signal_roi_freq = WireCell::Waveform::dft(signal_roi);
143  WireCell::Waveform::shrink(signal_roi_freq,respec);
144  for (size_t i=0;i!=signal_roi_freq.size();i++){
145  double freq;
146  // assuming 2 MHz digitization
147  if (i <signal_roi_freq.size()/2.){
148  freq = i/(1.*signal_roi_freq.size())*2.;
149  }else{
150  freq = (signal_roi_freq.size() - i)/(1.*signal_roi_freq.size())*2.;
151  }
152  std::complex<float> factor = filter_time(freq)*filter_low_loose(freq);
153  signal_roi_freq.at(i) = signal_roi_freq.at(i) * factor;
154  }
155  WireCell::Waveform::realseq_t signal_roi_decon = WireCell::Waveform::idft(signal_roi_freq);
156 
157  std::map<int, bool> flag_replace;
158  for (auto roi: rois){
159  flag_replace[roi.front()] = false;
160  }
161 
162  // judge if any ROI is good ...
163  for (auto roi: rois){
164  const int bin0 = std::max(roi.front()-1, 0);
165  const int binf = std::min(roi.back()+1, nbin-1);
166  double max_val = 0;
167  double min_val = 0;
168  for (int i=bin0; i<=binf; i++){
169  int time_bin = i-res_offset;
170  if (time_bin <0) time_bin += nbin;
171  if (time_bin >=nbin) time_bin -= nbin;
172  if (i==bin0){
173  max_val = signal_roi_decon.at(time_bin);
174  min_val = signal_roi_decon.at(time_bin);
175  }else{
176  if (signal_roi_decon.at(time_bin) > max_val) max_val = signal_roi_decon.at(time_bin);
177  if (signal_roi_decon.at(time_bin) < min_val) min_val = signal_roi_decon.at(time_bin);
178  }
179  }
180 
181  // if (signal.ch==1027)
182  //std::cout << roi.front() << " Xin " << max_val << " " << decon_limit1 << std::endl;
183 
184  if ( max_val > decon_limit1 && fabs(min_val) < max_val * roi_min_max_ratio)
185  flag_replace[roi.front()] = true;
186  }
187 
188  WireCell::Waveform::realseq_t temp_medians = medians;
189 
190  for (auto roi : rois) {
191  // original code used the bins just outside the ROI
192  const int bin0 = std::max(roi.front()-1, 0);
193  const int binf = std::min(roi.back()+1, nbin-1);
194  if (flag_replace[roi.front()]){
195  const double m0 = temp_medians[bin0];
196  const double mf = temp_medians[binf];
197  const double roi_run = binf - bin0;
198  const double roi_rise = mf - m0;
199  for (auto bin : roi) {
200  const double m = m0 + (bin - bin0)/roi_run*roi_rise;
201  temp_medians.at(bin) = m;
202  }
203  }
204  }
205 
206  // collection plane, directly subtracti ...
207  for (int i=0;i!=nbin;i++){
208  if (fabs(signal.at(i)) > 0.001) {
209  signal.at(i) = signal.at(i) - temp_medians.at(i) * scaling;
210  }
211  }
212 
213 
214  }else{
215  // collection plane, directly subtracti ...
216  for (int i=0;i!=nbin;i++){
217  if (fabs(signal.at(i)) > 0.001) {
218  signal.at(i) = signal.at(i) - medians.at(i) * scaling;
219  }
220  }
221  }
222  chansig[ch] = signal;
223  }
224 
225  // for (auto it: chansig){
226  // std::cout << "Xin2 " << it.second.at(0) << std::endl;
227  // }
228 
229  return true;
230 }
genie::units::m
static const double m
Definition: Units.h:79
WireCell::Waveform::realseq_t
Sequence< real_t > realseq_t
Definition: Waveform.h:31
WireCell::IChannelFilter::signal_t
Waveform::realseq_t signal_t
Definition: IChannelFilter.h:19
filter_low_loose
double filter_low_loose(double freq)
Definition: Microboone.cxx:50
WireCell::Waveform::sum2
Val sum2(const Sequence< Val > &seq)
Return sum of square of all entries in sequence.
Definition: Waveform.h:184
WireCell::SigProc::Derivations::CalcRMS
std::pair< double, double > CalcRMS(const WireCell::Waveform::realseq_t &signal)
Definition: Derivations.cxx:7
WireCell::Waveform::dft
compseq_t dft(realseq_t seq)
Definition: Waveform.cxx:141
WireCell::Waveform::shrink
void shrink(Sequence< Val > &seq, const Sequence< Val > &other)
Shrink (divide) seq values by values from the other sequence.
Definition: Waveform.h:162
max
static int max(int a, int b)
Definition: fortranscanner.cpp:60366
WireCell::Waveform::idft
realseq_t idft(compseq_t spec)
Definition: Waveform.cxx:149
cet::sqlite::min
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
bin
QTextStream & bin(QTextStream &s)
Definition: qtextstream.cpp:2006
WireCell::Waveform::compseq_t
Sequence< complex_t > compseq_t
A complex-valued sequence, eg for discrete spectrum powers.
Definition: Waveform.h:34
filter_time
double filter_time(double freq)
Definition: Microboone.cxx:35
Generated by   doxygen 1.8.11