RdChannelBeaconSignalExtractor.cc

Go to the documentation of this file.

#include "RdChannelBeaconSignalExtractor.h"
#include <fwk/CentralConfig.h>
#include <utl/ErrorLogger.h>

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

#include <utl/Trace.h>
#include <utl/TraceAlgorithm.h>
#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/config.h>
#include <utl/AugerUnits.h>
#include <utl/MathConstants.h>
#include <utl/TimeStamp.h>
#include <utl/TimeInterval.h>
#include <utl/FFTDataContainerAlgorithm.h>
#include <utl/AugerException.h>

#include <det/Detector.h>
#include <rdet/RDetector.h>

#include <complex>
#include <cmath>
#include <string>
#include <sstream>
#include <iostream>
#include <fstream>
#include <vector>
#include <map>

using namespace std;


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


namespace RdChannelBeaconSignalExtractor {

  RdChannelBeaconSignalExtractor::RdChannelBeaconSignalExtractor():
  fInfoLevel(eDebug),
  fBeaconFrequencies(vector<double>()),
  fReferenceChannel(1)
  {
  }


  RdChannelBeaconSignalExtractor::~RdChannelBeaconSignalExtractor()
  {
  }


  VModule::ResultFlag
  RdChannelBeaconSignalExtractor::Init()
  {
    // Initialize your module here. This method 
    // is called once at the beginning of the run.
    // The eSuccess flag indicates the method ended
    // successfully.  For other possible return types, 
    // see the VModule documentation.

    INFO("RdChannelBeaconSignalExtractor::Init()");

    // Read in the configurations of the xml file
    Branch topBranch =
      CentralConfig::GetInstance()->GetTopBranch("RdChannelBeaconSignalExtractor");
    topBranch.GetChild("infoLevel").GetData(fInfoLevel);
    if (fInfoLevel < eNone || fInfoLevel > eDebug) {
      ERROR("<infoLevel> out of bounds");
      return eFailure;
    }

    topBranch.GetChild("BeaconFrequencies").GetData(fBeaconFrequencies);
    topBranch.GetChild("ReferenceChannel").GetData(fReferenceChannel);

    return eSuccess;
  }


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

    // Check if there are events at all
    if(!event.HasREvent()) {
      WARNING("RdChannelBeaconSignalExtractor::No radio event found!");
      return eContinueLoop;
    } 
    
    // check if beacon frequencies have to be read from database, or are provided in XML file
    // in the first case, fBeaconFrequencies is empty, i.e. size()==0  
    if(fBeaconFrequencies.size() < 1) {
      fBeaconFrequencies = det::Detector::GetInstance().GetRDetector().GetBeaconFrequencies();
    }
    
    // Check if any frequencies are provided
    stringstream fMessage;
    fMessage.str("");
    fMessage << "Number of beacon frequencies to look at: "
             << fBeaconFrequencies.size();
    INFO(fMessage.str());

    if(fBeaconFrequencies.size() < 1) {
      WARNING("You must provide at least one frequency!");
      return eBreakLoop;
    }


    REvent& rEvent = event.GetREvent();

    // match time stamps of stations (if different)
    // this functionalty should be moved to its own module later
    matchStationTimeStamps(rEvent);

    double referenceAmplitude = 0; //  amplitude of first beacon frequency, of first station
    
    // loop through stations and for every station through every channel
    for (REvent::CandidateStationIterator sIt = rEvent.CandidateStationsBegin();
         sIt != rEvent.CandidateStationsEnd(); ++sIt) {
      Station& station = *sIt;
      
      for (Station::ChannelIterator cIt = station.ChannelsBegin();
           cIt != station.ChannelsEnd(); ++cIt) {
        Channel& channel = *cIt;
      
        // ignore all channels, but the reference channel (usually channel 1 = NS, high gain)
        if (channel.GetId()!=fReferenceChannel)
          continue;
              
        if (fInfoLevel>=eDebug) {              
          fMessage.str("");
          fMessage << "Reading data from station "
                    << channel.GetStationId()
                    << ", channel "
                    << channel.GetId()
                    << ".";
          INFO(fMessage.str());
        }  
            
        // Read the frequency spectrum of this channel
        revt::ChannelFrequencySpectrum& spectrum =
          channel.GetChannelFrequencySpectrum();

        set<int> beaconFrequencyBins; // can be different for each station because of different electronics and sampling frequencies
       
        // loop through beacon frequencies and find corresponding frequency bins in the spectrum 
        for (ChannelTimeSeries::SizeType i=0; i < fBeaconFrequencies.size(); ++i) {
          if (fInfoLevel>=eDebug) {              
            fMessage.str("");
            fMessage << "Beacon frequency: "
                    << fBeaconFrequencies[i]/megahertz
                    << " MHz";
            INFO(fMessage.str());

            fMessage.str("");
            fMessage << "Spectrum reaches from "
                    << channel.GetFrequencyOfBin(0)/megahertz
                    << " MHz to "
                    << channel.GetFrequencyOfBin(spectrum.GetSize()-1)/megahertz
                    << " MHz, binning = "
                    << spectrum.GetBinning()/megahertz
                    << " MHz.";
            INFO(fMessage.str());
          }  

          // check if spectrum has the correct lenght, such that the frequency bins are at the expected values
          // for which the beacon has been optimized
          if ( ((spectrum.GetSize()-1) % 1024) != 0 ) {
            fMessage.str("");
            fMessage << "Trace of station " << station.GetId() << " does not consist of 2048 samples or an integer multiple of it! "
                     << "It has " << (spectrum.GetSize()-1)*2 << " samples.";
            WARNING(fMessage.str());
            //return eFailure;
          }

          // check if frequency is inside the range (take care that the order of the spectrum is unknown)
          if( (fBeaconFrequencies[i] < min(channel.GetFrequencyOfBin(0),channel.GetFrequencyOfBin(spectrum.GetSize()-1)))
            || (fBeaconFrequencies[i] >  max(channel.GetFrequencyOfBin(0),channel.GetFrequencyOfBin(spectrum.GetSize()-1))) ) {
            WARNING("Beacon frequency must be inside of the spectrum!");
            return eContinueLoop;
          }

          // find nearest frequency bin:
          // bin = (f_beacon-f_0) * #lastBin / (f_1 - f_0)
          double bin = (fBeaconFrequencies[i]-channel.GetFrequencyOfBin(0)) * (spectrum.GetSize()-1)
                       / (channel.GetFrequencyOfBin(spectrum.GetSize()-1) - channel.GetFrequencyOfBin(0));
          beaconFrequencyBins.insert(round(bin));             

          // get the phase
          double phase = 0;

          phase = arg(spectrum[round(bin)]);
          if (fInfoLevel>=eDebug) {              
            fMessage.str("");
            fMessage << "Looking at frequency: "
                     << channel.GetFrequencyOfBin(round(bin)) / megahertz
                     << " MHz "
                     << "Phase: " << phase;
            INFO(fMessage.str());
          }
          
          if ((i==0) && (referenceAmplitude == 0))
            referenceAmplitude = abs(spectrum[round(bin)]);
          
        }
        
        if (fInfoLevel>=eDebug) {              
          fMessage.str("");
          fMessage << "Found: "
                   << beaconFrequencyBins.size()
                   << " beacon frequency bins.";
          INFO(fMessage.str());
        }
        if(beaconFrequencyBins.size() != fBeaconFrequencies.size()) {
          WARNING("Number of found beacon frequency bins differs from number of provided beacon frequencies.");
          return eBreakLoop;
        }
        
        // suppress spectrum outside of beacon frequencies,
        // and normalize amplitude of beacon frequencies to amplitude of first beacon frequency
        for (unsigned int i=0; i < spectrum.GetSize(); ++i) {
          if ( beaconFrequencyBins.find(i) == beaconFrequencyBins.end() ) {
            spectrum[i] *= 1.e-3;
          } else {
            spectrum[i] *= referenceAmplitude / abs(spectrum[i]);
            //cout << "Station " << channel.GetStationId() << " - Spectrum value at beacon freq # " << i+1 << " = " << spectrum[i] << endl;
          }          
        }
      }
    }

    return eSuccess;
  }


  VModule::ResultFlag
  RdChannelBeaconSignalExtractor::Finish()
  {
    stringstream fMessage;
    // Put any termination or cleanup code here.
    // This method is called once at the end of the run.

    INFO("RdChannelBeaconSignalExtractor::Finish()");
 
    return eSuccess;
  }

  void
  RdChannelBeaconSignalExtractor::matchStationTimeStamps(revt::REvent& rEvent) const
  {
    // store time stamp of the first station as reference
    double firstTimeStamp  = rEvent.CandidateStationsBegin()->GetRecData().GetParameter(eTraceStartTime);
    // loop through all stations:
    // shift traces if the time stamp does not match the one of the first station
    for (REvent::CandidateStationIterator sIt = rEvent.CandidateStationsBegin();
         sIt != rEvent.CandidateStationsEnd(); ++sIt) {
      Station& station = *sIt;
      double stationTimeStamp = station.GetRecData().GetParameter(eTraceStartTime);

      if (stationTimeStamp != firstTimeStamp) {
        // correct for the difference and set new time stamp
        for (Station::ChannelIterator cIt = station.ChannelsBegin();
             cIt != station.ChannelsEnd(); ++cIt){
            Channel& channel = *cIt;
        
            FFTDataContainerAlgorithm::ShiftTimeSeries(channel.GetFFTDataContainer(), stationTimeStamp-firstTimeStamp);
        }
        station.GetRecData().SetParameter(eTraceStartTime,firstTimeStamp);
      }
    }
  }

}