RdDummyChannelToStationConverter.cc

Go to the documentation of this file.

#include "RdDummyChannelToStationConverter.h"

#include <fwk/CentralConfig.h>

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

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



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

namespace RdDummyChannelToStationConverter {

  RdDummyChannelToStationConverter::RdDummyChannelToStationConverter()
  {
  }


  RdDummyChannelToStationConverter::~RdDummyChannelToStationConverter()
  {
  }


  VModule::ResultFlag
  RdDummyChannelToStationConverter::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("RdDummyChannelToStationConverter::Init()");
    Branch topBranch =
      CentralConfig::GetInstance()->GetTopBranch("RdDummyChannelToStationConverter");

    topBranch.GetChild("UseIntegratedAntennaPattern").GetData(fUseIntegratedAntennaPattern);

    return eSuccess;
  }


  VModule::ResultFlag
  RdDummyChannelToStationConverter::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();
//  stringstream fMessage;
    Vector3D VDZero;
    VDZero = 0.0, 0.0, 0.0;
    stringstream fMessage;
    // loop through stations and for every station through every channel
    for (REvent::StationIterator sIt = rEvent.StationsBegin();
         sIt != rEvent.StationsEnd(); ++sIt) {
      Station& station = *sIt;
      StationTimeSeries& newstationtimeseries=station.GetStationTimeSeries();
      bool firststationchannel = true;

      for (Station::ChannelIterator cIt = station.ChannelsBegin();
           cIt != station.ChannelsEnd(); ++cIt) {
        Channel& channel = *cIt;
        //only works for active channels
        if (!channel.IsActive())
          continue;

        fMessage.str("");
        fMessage << "Copying channel "
               << channel.GetId()
               << " of station "
               << channel.GetStationId()
               << " to the StationTimeSeries.";
        INFO(fMessage.str());

        // Work on the trace of this channel
        const revt::ChannelTimeSeries& channeltrace =
          channel.GetChannelTimeSeries();

        // when working on the first channel, fill station timeseries with zeroes and set some station-specific quantities
        if (firststationchannel) {
          firststationchannel = false;
          // overwrite the station time series with an empty one of the correct size
          newstationtimeseries = StationTimeSeries(channeltrace.GetSize(),channeltrace.GetBinning(),VDZero);
          // set the Nyquist zone of the station, does no check if Nyquist zones of Channels are the same!
          station.SetNyquistZone(channel.GetNyquistZone());
        }

        // check if newstationtimeseries.GetSize() is correct
        if (newstationtimeseries.GetSize() != channeltrace.GetSize()) {
          fMessage.str("");
          fMessage << "Number of samples in trace of "
                   << channel.GetId()
                   << " belonging to station "
                   << channel.GetStationId()
                   << " is different than the number of samples in other channels of that station.";
          WARNING(fMessage.str());
          return eFailure;
        }

        const rdet::RDetector& rDetector =
          det::Detector::GetInstance().GetRDetector();
        const rdet::Channel& detChannel =
          rDetector.GetStation(station.GetId()).GetChannel(channel.GetId());
        double orientAzimuth = detChannel.GetOrientationAzimuth();
        double orientZenith = detChannel.GetOrientationZenith();
        
        if (orientAzimuth < 0 && orientAzimuth > -360 * deg)
          orientAzimuth += 360 * deg;

        // Doing some kind of dodgy association of the channels to the coordinates
        // The design-decision is that this module is only for testing purposes so
        // a proper coordinate transformation is not done here.
        fMessage.str("");
        double alignMargin = 5*deg;
        if ( orientZenith > alignMargin ) {
          fMessage << "Encountering a channel that seems to be vertically aligned, not converting.";
          WARNING(fMessage.str());
          continue;
        }

        double correspondingCoordinate; // Will always be initialized or ignored in the code below
        double correspondingSign; // Will always be initialized or ignored in the code below
        fMessage.str("");
        fMessage << "Associating channel " << channel.GetId() << " for station " << station.GetId() << " at " <<
          orientAzimuth/deg << " deg azimuthal angle";
        if ( orientAzimuth < alignMargin or orientAzimuth > 360 * deg - alignMargin ) {
          fMessage << " with a positive x coordinate.";
          INFO(fMessage.str());
          correspondingCoordinate = 0;
          correspondingSign = 1;
        } else if ( orientAzimuth > 90 * deg - alignMargin or orientAzimuth < 90 * deg + alignMargin ) {
          fMessage << " with a positive y coordinate.";
          INFO(fMessage.str());
          correspondingCoordinate = 1;
          correspondingSign = 1;
        } else if ( orientAzimuth > 180 * deg - alignMargin or orientAzimuth < 180 * deg + alignMargin ) {
          fMessage << " with a negative x coordinate.";
          INFO(fMessage.str());
          correspondingCoordinate = 0;
          correspondingSign = -1;
        } else if ( orientAzimuth > 270 * deg - alignMargin or orientAzimuth < 270 * deg + alignMargin ) {
          fMessage << " with a negative y coordinate.";
          INFO(fMessage.str());
          correspondingCoordinate = 1;
          correspondingSign = -1;
        } else {
          fMessage.str("");
          fMessage << "No good dodgy conversion possible because the antenna does not seem to be aligned with " <<
            "any spatial coordinate within " << alignMargin/deg << "deg.";
          WARNING(fMessage.str());
          continue;
        }

        // if UseIntegratedAntennaResponse is active, scale the Channel spectrum with the integrated antenna response
        if (fUseIntegratedAntennaPattern) {
          // retrieve a write-access reference of the channel spectrum
          revt::ChannelFrequencySpectrum& channelspectrum = channel.GetChannelFrequencySpectrum();
          // get the design bandwidth from the channel
          const double lowerDesignFreq = detChannel.GetDesignLowerFreq();
          const double upperDesignFreq = detChannel.GetDesignUpperFreq();

          // get the integrated antenna response and apply it to the data
          for (ChannelFrequencySpectrum::SizeType i=0; i<channelspectrum.GetSize(); ++i) {
            const double freq = channel.GetFrequencyOfBin(i);
            // apply antenna response, if outside design band set to zero
            if ((freq >= lowerDesignFreq) && (freq <= upperDesignFreq))
              channelspectrum[i] /= detChannel.GetIntegratedEffectiveAntennaHeight(freq);
            else
              channelspectrum[i] = 0.0;
          }
        }

        // if UseIntegratedAntennaResponse is active, the frequency spectrum has been changed
        // we therefore have to make sure the time domain representation is up-to-date
        // consequently, we request a new reference to the channel trace
        const revt::ChannelTimeSeries& newchanneltrace = channel.GetChannelTimeSeries();
        for (ChannelTimeSeries::SizeType i=0; i<newchanneltrace.GetSize(); ++i) {
          newstationtimeseries[i][correspondingCoordinate] = correspondingSign * newchanneltrace[i];
        }
      }
    }



    return eSuccess;
  }


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

    INFO("RdDummyChannelToStationConverter::Finish()");

    return eSuccess;
  }

}