RdChannelBeaconSimulator.cc

Go to the documentation of this file.

#include "RdChannelBeaconSimulator.h"

#include <fwk/CentralConfig.h>
#include <fwk/RunController.h>
#include <fwk/LocalCoordinateSystem.h>
#include <fwk/CoordinateSystemRegistry.h>
#include <fwk/VModule.h>

#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/config.h>
#include <utl/TraceAlgorithm.h>
#include <utl/PhysicalConstants.h>
#include <utl/AxialVector.h>
#include <utl/Math.h>
#include <utl/FFTDataContainerAlgorithm.h>

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

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

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

namespace RdChannelBeaconSimulator {

  template<typename G>
  inline
  string
  ToString(const G& g, const CoordinateSystemPtr cs)
  {
    ostringstream os;
    os << '('
       << g.GetX(cs)/meter << ", "
       << g.GetY(cs)/meter << ", "
       << g.GetZ(cs)/meter << ") [m]";
    return os.str();
  }

  RdChannelBeaconSimulator::RdChannelBeaconSimulator():
  fBeaconFrequencies(vector<double>())
  {
  }

  RdChannelBeaconSimulator::~RdChannelBeaconSimulator()
  {
  }

  VModule::ResultFlag
  RdChannelBeaconSimulator::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("RdChannelBeaconSimulator::Init()");

    // Read in the configurations of the xml file
    Branch topBranch =
      CentralConfig::GetInstance()->GetTopBranch("RdChannelBeaconSimulator");
    topBranch.GetChild("BeaconFrequencies").GetData(fBeaconFrequencies);
    topBranch.GetChild("BeaconPosition").GetData(fBeaconPosition);
    topBranch.GetChild("Amplitude").GetData(fAmplitude);
    topBranch.GetChild("TimeJitter").GetData(fTimeJitter);
    topBranch.GetChild("TimeJitterMean").GetData(fTimeJitterMean);
    topBranch.GetChild("OutputShiftFile").GetData(fOutputShiftFile);

    return eSuccess;
  }

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

    fMessage.str("");
    fMessage << "Start Beaconsimulator!!! ";
    INFO(fMessage.str());

    REvent& rEvent = event.GetREvent();
    // get detector description (for antenna positions)
    const det::Detector& detector = det::Detector::GetInstance();
    const rdet::RDetector& rDetector = detector.GetRDetector();
    // get reference coordinate system of detector (usually PampaAmarilla)
    const CoordinateSystemPtr referenceCS = detector.GetReferenceCoordinateSystem();

    // initialize random number generator
    srand( time(NULL));

    float shift = 0.0;
    ofstream outputfile;

    if(fOutputShiftFile != "") {
      outputfile.open(fOutputShiftFile.c_str(), ofstream::app);
    }

    // the absolute points in time and space
    const Point siteOrigin(0,0,0, referenceCS);
    Vector vStationPosition(0,0,0, referenceCS);
    Vector vBeaconPosition(fBeaconPosition[0],fBeaconPosition[1],fBeaconPosition[2], referenceCS);
    Vector vPositionDifference(0,0,0, referenceCS);

    sleep(1);
    // loop through stations and for every station through every channel
    for (REvent::StationIterator sIt = rEvent.StationsBegin();
         sIt != rEvent.StationsEnd(); ++sIt) {
      Station& station = *sIt;

      //First try to get the position
      const rdet::Station& dStation = rDetector.GetStation(*sIt);
      vStationPosition = dStation.GetPosition() - siteOrigin;

      //Calculate the vector and the distance from Beacon to the station
      vPositionDifference = vBeaconPosition - vStationPosition;
      double distance = sqrt(vPositionDifference * vPositionDifference);

      fMessage.str("");
      fMessage << "Position of Station " << dStation.GetId() << ": " << ToString(vStationPosition,  referenceCS)<< endl
               << "Difference to Beacon: " << ToString(vPositionDifference, referenceCS)
               << " => Distance to Beacon: " << distance << " [m]";
      INFO(fMessage.str());

      if (fTimeJitter == 1) {
        // randomly shift the timeseries to emulate the GPS jitter
//        float shift = ((rand() / (float)RAND_MAX) * 6. ) - 3.;
        // try gaussian distibution
        shift = (1/sqrt(2)) * fTimeJitterMean * sqrt( -2. * log((float)rand() / ((float)RAND_MAX + 1.))) * sin( 2. * kPi * ((float)rand() / ((float)RAND_MAX + 1.)) );
      }
      // loop through channels
      for (Station::ChannelIterator cIt = station.ChannelsBegin();
           cIt != station.ChannelsEnd(); ++cIt) {
        Channel& channel = *cIt;
        if (!channel.IsActive())
          continue;

        // Read the frequency spectrum of this channel
        revt::ChannelFrequencySpectrum& spectrum =
          channel.GetChannelFrequencySpectrum();

        // loop through beacon frequencies and find corresponding frequency in the spectrum
        for (ChannelTimeSeries::SizeType i=0; i < fBeaconFrequencies.size(); ++i) {

          //Calculate the distance in wavelength and the corresponding phase at the station
          double n_wav = ((fBeaconFrequencies[i]/hertz)*distance)/300000000. ;
          int result = static_cast<int>(n_wav/1.);
          double mod = n_wav - static_cast<double>(result);
          double phase = mod * kTwoPi;
//          double amplitude = 10000. * pow(distance/1000, -4);
          double amplitude = fAmplitude * pow(distance/1000, -2);

          fMessage.str("");
          fMessage << "Beacon frequency: " << fBeaconFrequencies[i]/megahertz << " MHz"
                   << " Phase: " << phase << " Amplitude: " << amplitude;
          INFO(fMessage.str());

          // 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));

          //Sum up calculated amplitude and original amplitude and set bin
          spectrum[bin] += amplitude * exp(complex<double>(0,1)*phase);

          // get the phase and amplitude of the bin
          phase = 0;
          amplitude = 0;

          phase = arg(spectrum[round(bin)]);
          amplitude = abs(spectrum[round(bin)]);
        }

        if(fOutputShiftFile != "") {
          if ( channel.GetId() == 1 )
          {
            outputfile << shift << endl;
          }
        }
        //Shift the trace of the channel
        FFTDataContainerAlgorithm::ShiftTimeSeries(channel.GetFFTDataContainer(), shift);
      }
    }
    if(fOutputShiftFile != "")
      outputfile.close();
    return eSuccess;
  }

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

    INFO("RdChannelBeaconSimulator::Finish()");

    return eSuccess;
  }

}