RdChannelSineWaveSuppressor.cc

Go to the documentation of this file.

#include "RdChannelSineWaveSuppressor.h"
#include <utl/ErrorLogger.h>

#include <fwk/CentralConfig.h>

#include <evt/Event.h>
#include <revt/REvent.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerRRecData.h>
#include <revt/Header.h>
#include <revt/Station.h>
#include <revt/Channel.h>
#include <revt/StationRecData.h>

#include <rdet/RDetector.h>
#include <rdet/Channel.h>

#include <utl/Trace.h>
#include <utl/TraceAlgorithm.h>
#include <utl/Reader.h>
#include <utl/config.h>
#include <utl/AugerUnits.h>
#include <utl/Math.h>

#include <TMath.h>
#include <TMinuit.h>

using namespace std;
using namespace utl;
using namespace fwk;
using namespace revt;
using namespace rdet;


namespace RdChannelSineWaveSuppressor {

  REvent* RdChannelSineWaveSuppressor::fgCurrentREvent;
  int RdChannelSineWaveSuppressor::fgCurrentStation;
  int RdChannelSineWaveSuppressor::fgCurrentChannel;
  int RdChannelSineWaveSuppressor::fgFitWindow;
  unsigned int RdChannelSineWaveSuppressor::fgSignalStart;
  unsigned int RdChannelSineWaveSuppressor::fgSignalStop;
  unsigned int RdChannelSineWaveSuppressor::fgNoiseStart;
  unsigned int RdChannelSineWaveSuppressor::fgNoiseStop;


  VModule::ResultFlag
  RdChannelSineWaveSuppressor::Init()
  {
    INFO("RdChannelSineWaveSuppressor::Init()");

    // Read in the configurations of the xml file
    Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdChannelSineWaveSuppressor");
    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);

    topBranch.GetChild("ScanWindowSize").GetData(fScanWindowSize);
    topBranch.GetChild("NoiseInterval").GetData(fNoiseInterval);
    topBranch.GetChild("FitWindow").GetData(fgFitWindow);

    if (fNoiseInterval >= 1.0) {
      ERROR("<NoiseInterval> must be a fraction smaller than one");
      return eFailure;
    }

    if (fgFitWindow > 3) {
      ERROR("<FitWindow> out of bounds");
      return eFailure;
    }
    return eSuccess;
  }


  VModule::ResultFlag
  RdChannelSineWaveSuppressor::Run(evt::Event& event)
  {
    INFO("RdChannelSineWaveSuppressor::Run()");

    REvent& rEvent = event.GetREvent();
    fgCurrentREvent = &rEvent;
    ostringstream info;

    // Loop over all stations and channels
    for (auto& station : rEvent.StationsRange()) {

      fgCurrentStation = station.GetId();

      for (auto& channel : station.ChannelsRange()) {

        if (!channel.IsActive())
          continue;

        fgCurrentChannel = channel.GetId();

        const revt::ChannelFrequencySpectrum& spectrum = channel.GetChannelFrequencySpectrum();
        double spectrumBinning = spectrum.GetBinning();
        const rdet::RDetector& rDetector = det::Detector::GetInstance().GetRDetector();
        const rdet::Channel& detChannel =
          rDetector.GetStation(station.GetId()).GetChannel(channel.GetId());
        const double lowerDesignFreq = detChannel.GetDesignLowerFreq();
        const double upperDesignFreq = detChannel.GetDesignUpperFreq();
        unsigned int lowerBin = (int)(lowerDesignFreq/spectrumBinning);
        unsigned int upperBin = (int)(upperDesignFreq/spectrumBinning);
        int WindowSize;

        if (fScanWindowSize > upperBin-lowerBin) {
          ostringstream warn;
          warn << "<ScanWindowSize> too big, window size set to full spectrum width ("
               << (upperBin - lowerBin) << " bins)";
          WARNING(warn);
          WindowSize = upperBin - lowerBin;
        } else
          WindowSize = fScanWindowSize;

        const unsigned int centralBin = WindowSize/2 + 1;
        vector<double> vSlidingWindow;
        vector<double> vThreshold;
        vector<double> vMaxLikelihood;
        vector<int> vRFIBins;
        double RFIThreshold = 0;
        double squaredSum = 0;
        double maxLikelihood = 0;

        // Caluclate argument of CDF on sample level assuming Rayleigh distribution of noise spectrum
        double InvSampleProp = 1 - pow(fNoiseInterval, 1. / (upperBin - lowerBin));
        double CDFarg = -log(InvSampleProp);

        // Move sliding window over frequency spectrum to find RFI frequencies
        for (ChannelFrequencySpectrum::SizeType i = lowerBin; i < lowerBin+WindowSize; ++i) {
          vSlidingWindow.push_back(abs(spectrum[i]));
          squaredSum += Sqr(abs(spectrum[i]));
        }

        maxLikelihood = sqrt(squaredSum / (2*WindowSize));
        RFIThreshold = sqrt(CDFarg * 2 * Sqr(maxLikelihood));

        for (vector<double>::size_type i = 0; i <= centralBin; ++i) {
          if (vSlidingWindow[i] > RFIThreshold) {
            vRFIBins.push_back(lowerBin + i);
            vThreshold.push_back(RFIThreshold);
            vMaxLikelihood.push_back(maxLikelihood);
          }
        }

        for (ChannelFrequencySpectrum::SizeType i = lowerBin+WindowSize; i <= upperBin; ++i) {
          squaredSum += Sqr(vSlidingWindow.front()) + Sqr(abs(spectrum[i]));
          maxLikelihood = sqrt(squaredSum / (2 * WindowSize));
          RFIThreshold = sqrt(CDFarg * 2 * Sqr(maxLikelihood));
          vSlidingWindow.erase(vSlidingWindow.begin());
          vSlidingWindow.push_back(abs(spectrum[i]));

          if (vSlidingWindow[centralBin] > RFIThreshold) {
            vRFIBins.push_back(i - WindowSize + centralBin + 1);
            vThreshold.push_back(RFIThreshold);
            vMaxLikelihood.push_back(maxLikelihood);
          }
        }

        for (vector<double>::size_type i = centralBin; i < vSlidingWindow.size(); ++i) {
          if (vSlidingWindow[i] > RFIThreshold) {
            vRFIBins.push_back(upperBin - WindowSize + i + 1);
            vThreshold.push_back(RFIThreshold);
            vMaxLikelihood.push_back(maxLikelihood);
          }
        }

        // Just suppose nothing is above threshold, continue with next channel
        if (vRFIBins.empty()) {
          info << "No RFI identified in channel " << channel.GetId();
          INFOIntermediate(info);
          continue;
        }

        int PeakWidth = 1;
        int PeakBin;
        vector<int> vFitBins;

        // Select highest bins from RFI sources wider than one bin
        for (vector<int>::size_type i = 0; i < vRFIBins.size(); ++i) {
          if (i < vRFIBins.size() - 1) {
            if (vRFIBins[i+1] - vRFIBins[i] == 1)
              ++PeakWidth;
            else {
              PeakBin = vRFIBins[i - PeakWidth + 1];
              for (vector<int>::size_type j = i - PeakWidth + 1; j <= i; ++j) {
                if (abs(spectrum[vRFIBins[j]]) > abs(spectrum[PeakBin]))
                  PeakBin = vRFIBins[j];
              }
              vFitBins.push_back(PeakBin);
              PeakWidth = 1;
            }
          } else { // Take care of last element
            PeakBin = vRFIBins[i - PeakWidth + 1];
            for (vector<int>::size_type j = i - PeakWidth + 1; j <= i; ++j) {
              if (abs(spectrum[vRFIBins[j]]) > abs(spectrum[PeakBin]))
                  PeakBin = vRFIBins[j];
            }
            vFitBins.push_back(PeakBin);
          }
        }

        // Sort RFI sources to fit strongest contribution first
        if (vFitBins.size() > 1) {
          RdChannelSineWaveSuppressor::RFISort(vFitBins, 0, vFitBins.size() - 1);
        }

        info << "Identified RFI frequencies in channel " << channel.GetId()
             << " of station " << station.GetId() << " in decreasing order of strength:";
        INFOIntermediate(info);

        for (auto& elem : vFitBins) {
          info << elem * spectrumBinning/MHz << " MHz ";
          INFODebug(info);
        }
        INFODebug("\n");

        revt::ChannelTimeSeries& trace = channel.GetChannelTimeSeries();
        double samplingFrequency = 1. / trace.GetBinning();

        // Set fitting globals defining windows
        fgSignalStart = station.GetRecData().GetParameter(revt::eSignalSearchWindowStart)*samplingFrequency;
        fgSignalStop = station.GetRecData().GetParameter(revt::eSignalSearchWindowStop)*samplingFrequency;
        fgNoiseStart = station.GetRecData().GetParameter(revt::eNoiseWindowStart)*samplingFrequency;
        fgNoiseStop = station.GetRecData().GetParameter(revt::eNoiseWindowStop)*samplingFrequency;

        // Loop over all RFI frequencies and perform the fit for each of them
        for (unsigned int j = 0; j < vFitBins.size(); ++j) {

          // Set fit threshold for fitting fequency and find noise level in vicinity of RFI bin
          double FitThreshold = 0;
          double FitNoise = 0;
          for (vector<int>::size_type i = 0; i < vRFIBins.size(); ++i) {
            if (vFitBins[j] == vRFIBins[i]) {
              FitThreshold = vThreshold[i];
              FitNoise = vMaxLikelihood[i];
              break;
            }
          }

          // Calculate fit limits and estimate fit values
          double FitFreq = vFitBins[j] * spectrumBinning;
          double Freqmin = FitFreq - 0.5*spectrumBinning;
          double Freqmax = FitFreq + 0.5*spectrumBinning;

          // spectrum corrected for noise level, estimate when frequency is in middle of bin, so mimimum leakage
          double AmpEst = (abs(spectrum[vFitBins[j]]) - FitNoise)*spectrumBinning*2.0;
          // average of noise floor (based on Rayleigh distribution)
          double Ampmin = FitNoise * spectrumBinning;
          //spectrum corrected for noise level, maximum spectral leakage
          double Ampmax = abs(spectrum[vFitBins[j]]) * spectrumBinning*TMath::Pi() - FitNoise*spectrumBinning*2.0;

          double PhaseEst = arg(spectrum[vFitBins[j]]) + FitFreq/samplingFrequency*2*TMath::Pi();

          TMinuit minuit(3);
          double argList[10];
          int errFlag = 0;
          argList[0] = -1;

          if (fInfoLevel >= eInfoDebug) {
            minuit.mnexcm("SET PRINTOUT", argList, 1, errFlag);
            minuit.mnexcm("SET NOWARNINGS", argList, 0, errFlag);
          }

          // Set minimization function and parameters
          minuit.SetFCN(RdChannelSineWaveSuppressor::SineFitFnc);
          argList[0] = 1;
          minuit.mnexcm("SET ERRORDEF", argList, 1, errFlag);
          minuit.mnparm(0, "amp", AmpEst, 0.00001, Ampmin, Ampmax, errFlag);
          minuit.mnparm(1, "freq", FitFreq, 0.000001, Freqmin, Freqmax, errFlag);
          minuit.mnparm(2, "phase", PhaseEst, 0.00001, 0, 0, errFlag);

          // Perform minimization
          argList[0] = 5;
          minuit.mnexcm("CALI", argList, 1, errFlag);
          argList[0] = 500;
          minuit.mnexcm("MINIMIZE", argList, 1, errFlag);

          // Get fit results (errors not implemented yet)
          double amp = 0;
          double freq = 0;
          double phase = 0;
          double foo = 0;

          minuit.GetParameter(0, amp, foo);
          minuit.GetParameter(1, freq, foo);
          minuit.GetParameter(2, phase, foo);

          // Substract fitted wave from time series if fitted amplitude large enough
          if (amp > FitThreshold * spectrumBinning) {
            for (ChannelTimeSeries::SizeType i = 0; i < trace.GetSize(); ++i) {
              trace[i] = trace[i] - amp*sin(i/samplingFrequency*freq*2*TMath::Pi() + phase);
            }
            info << "Substracted sine wave from channel " << channel.GetId()
                 << " of station " << station.GetId() << " with frequency "
                 << freq/MHz << " MHz, amplitude " << amp*1e6
                 << " muV/m and phase " << phase/TMath::Pi() << " pi.";
            INFODebug(info);
          } else {
            info << "Fit failed for RFI source from channel " << channel.GetId()
                 << " of station " << station.GetId() << " with frequency "
                 << FitFreq/MHz << " MHz. Fitted amplitude of " << amp*1e6
                 << " muV/m is lower than RFI threshold of "
                 << FitThreshold*spectrumBinning*1e6 << " muV/m.";
            INFODebug(info);
          }
        }
      }
    }
    return eSuccess;
  }


  VModule::ResultFlag
  RdChannelSineWaveSuppressor::Finish()
  {
    INFO("RdChannelSineWaveSuppressor::Finish()");

    return eSuccess;
  }


  void
  RdChannelSineWaveSuppressor::SineFitFnc(int& /*nPar*/, double* const /*grad*/, double& value,
                                          double* const par, const int /*flag*/)
  {
    const double& amp = par[0];
    const double& freq = par[1];
    const double& phase = par[2];
    double Chi2 = 0;
    const revt::ChannelTimeSeries& trace =
        fgCurrentREvent->GetStation(fgCurrentStation).GetChannel(fgCurrentChannel).GetChannelTimeSeries();
    double samplingFrequency = 1./trace.GetBinning();

    for (ChannelTimeSeries::SizeType i = 0; i < trace.GetSize(); ++i) {
      if ((i > fgSignalStart && i < fgSignalStop && fgFitWindow > 1) ||
          (i > fgNoiseStart && i < fgNoiseStop && fgFitWindow == 3) ||
          (i < fgNoiseStart && i > fgNoiseStop && fgFitWindow == 1))
        continue;

      Chi2 += (trace[i] - amp*sin(i/samplingFrequency*freq*2*TMath::Pi() + phase)) *
              (trace[i] - amp*sin(i/samplingFrequency*freq*2*TMath::Pi() + phase));
    }
    value = Chi2;
  }

  void
  RdChannelSineWaveSuppressor::RFISort(std::vector<int> &v, int left, int right)
  {
    const revt::ChannelFrequencySpectrum& spectrum =
      fgCurrentREvent->GetStation(fgCurrentStation).GetChannel(fgCurrentChannel).GetChannelFrequencySpectrum();
    int i = left, j = right;
    int tmp;
    double pivot = abs(spectrum[v[(left + right) / 2]]);

    while (i <= j) {
      while (abs(spectrum[v[i]]) > pivot)
        ++i;
      while (abs(spectrum[v[j]]) < pivot)
        --j;
      if (i <= j) {
        tmp = v[i];
        v[i] = v[j];
        v[j] = tmp;
        ++i;
        --j;
      }
    }
    if (left < j)
      RdChannelSineWaveSuppressor::RFISort(v, left, j);
    if (i < right)
      RdChannelSineWaveSuppressor::RFISort(v, i, right);
  }

}