RdADCTraceFixer.cc

Go to the documentation of this file.

#include "RdADCTraceFixer.h"

#include <fwk/CentralConfig.h>

#include <utl/config.h>
#include <utl/ErrorLogger.h>
#include <utl/Reader.h>

#include <evt/Event.h>
#include <revt/REvent.h>
#include <revt/Station.h>
#include <revt/Channel.h>


using namespace revt;
using namespace fwk;
using namespace utl;
using std::ostringstream;


namespace RdADCTraceFixer {


  VModule::ResultFlag
  RdADCTraceFixer::Init()
  {
    Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdADCTraceFixer");
    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
    topBranch.GetChild("CorrectChannelMixing").GetData(fCorrectChannelMixing);

    return eSuccess;
  }


  VModule::ResultFlag
  RdADCTraceFixer::Run(evt::Event& event)
  {
    // Check if there are events at all
    if (!event.HasREvent()) {
      WARNING("No radio event found!");
      return eContinueLoop;
    }

    REvent& rEvent = event.GetREvent();

    // loop through stations and for every station through every channel
    for (auto& station : rEvent.StationsRange()) {

      //correct channel mixing (flaws in EA data from 11/22-03/23)
      if (fCorrectChannelMixing) {
        std::vector<std::pair<ChannelADCTimeSeries, ChannelADCTimeSeries>> traceOpts;
        std::pair<ChannelADCTimeSeries, ChannelADCTimeSeries> currTraces;

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

          if (!channel.IsActive())
            break;

          if (channel.GetId() == 0) {
            currTraces.first = channel.GetChannelADCTimeSeries();
          } else {
            currTraces.second = channel.GetChannelADCTimeSeries();
          }
        }

        const double binning = station.GetChannel(0).GetChannelADCTimeSeries().GetBinning();
        const int traceLength = currTraces.first.GetSize();

        // we test 3 options for the NS channel (one channel is conclusive):
        // no channel mixing
        traceOpts.push_back(currTraces);

        // trace starts on odd bin
        for (int i = 0; i < traceLength; i += 2) {
          currTraces.first[i] = traceOpts[0].first[i];
          currTraces.first[i+1] = traceOpts[0].second[i];
          currTraces.second[i] = traceOpts[0].second[i+1];
          currTraces.second[i+1] = traceOpts[0].first[i+1];
        }
        traceOpts.push_back(currTraces);

        // trace starts on even bin
        for (int i = 0; i < traceLength; i += 2) {
          currTraces.first[i] = traceOpts[0].first[i+1];
          currTraces.first[i+1] = traceOpts[0].second[i+1];
          currTraces.second[i] = traceOpts[0].second[i];
          currTraces.second[i+1] = traceOpts[0].first[i];
        }
        traceOpts.push_back(currTraces);

        // now test all options by making a FFT
        // look at frequency spectrum to decide which option is correct
        std::vector<double> freqRatios;
        for (auto const& opt : traceOpts) {
          Channel& channelNS = station.GetChannel(0);
          // we change timeTrace, which changes GetChannelTimeSeries(), which will update
          // the FFT when calling GetFrequencySpectrum()
          auto& timeTrace = channelNS.GetChannelTimeSeries();

          timeTrace.Clear();
          timeTrace.SetBinning(binning);
          timeTrace = opt.first * 1.0;

          auto& myFFTDataContainer = channelNS.GetFFTDataContainer();
          const auto& freqSpectrum = myFFTDataContainer.GetFrequencySpectrum();

          const double ratio = CalculateRatio(freqSpectrum, traceLength);
          freqRatios.push_back(ratio);
        }

        const int maxIndex = std::max_element(freqRatios.begin(), freqRatios.end()) - freqRatios.begin();

        ostringstream info;
        info <<"\n" << (maxIndex ? "Original trace correct " : "Channel mixing ")
             << "identified. Saved correct trace.";
        INFO(info);

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

          if (!channel.IsActive())
            break;

          auto& adcTrace = channel.GetChannelADCTimeSeries();

          if (channel.GetId() == 0) {
            adcTrace = traceOpts[maxIndex].first;
          } else {
            adcTrace = traceOpts[maxIndex].second;
          }
        }
      }
    }

    return eSuccess;
  }


  VModule::ResultFlag
  RdADCTraceFixer::Finish()
  {
    return eSuccess;
  }


  double
  RdADCTraceFixer::CalculateRatio(const ChannelFrequencySpectrum& freqSpectrum, int traceLength)
  {
    double freq = 0;
    double inRange = 0;
    double outRange = 0;
    int counter = 1;
    for (const auto& amp : freqSpectrum) {
      freq = double(counter) / double(traceLength) * 250;

      if (freq >= 30 && freq <= 80) {
        inRange += abs(amp);
      } else {
        outRange += abs(amp);
      }
      counter++;
    }
    return inRange / outRange;
  }

}