RdStationInterpolator.cc

Go to the documentation of this file.

#include "RdStationInterpolator.h"
#include <fwk/CentralConfig.h>
#include <fwk/VModule.h>

#include <rdet/RDetector.h>
#include <rdet/Station.h>
#include <utl/CoordinateSystem.h>
#include <utl/Point.h>
#include <evt/RadioSimulation.h>
#include <utl/TimeStamp.h>
#include <utl/TimeInterval.h>
#include <utl/UTCDateTime.h>
#include <utl/Trace.h>
#include <utl/AugerUnits.h>
#include <evt/Event.h>
#include <evt/ShowerSimData.h>
#include <revt/REvent.h>
#include <revt/EventTrigger.h>
#include <revt/Station.h>
#include <revt/Header.h>
#include <revt/StationTriggerData.h>
#include <revt/StationSimData.h>
#include <utl/StringCompare.h>
#include <cmath>
#include <algorithm>

#include <utl/FFTDataContainer.h>
#include <string>
#include <evt/SimRadioPulse.h>
#include <utl/MathConstants.h>
#include <complex>
#include <map>
#include <vector>
#include <iomanip>
#include <fstream>
#include <math.h>
#include <boost/utility.hpp>
#include <TRandom2.h>
#include <boost/algorithm/string.hpp>
#include <stdlib.h>
#include <rdet/Channel.h>


using namespace std;
using namespace evt;
using namespace det;
using namespace utl;
using namespace revt;
using namespace fwk;

#define OUT(x) if ((x) <= fInfoLevel) cerr << "  "


namespace RdStationInterpolator {

  typedef multimap<double, const SimRadioPulse&> PulseMMap;
  typedef pair<double, const SimRadioPulse&> PulsePair;


  RdStationInterpolator::RdStationInterpolator() :
    fMaximumDistance(5.*meter), // not yet used since not jet tested, up to which distance the interpolation works properly
    fTracePaddingAtFront(2000.*nanosecond),
    fMinimumTraceLength(12000.*nanosecond),
    fEventRate(0.0001*hertz),
    fEventTimeOffset(0.0),
    fUseRateTiming(false),
    //fIncludedStationIds,
    //fExcludedStationIds
    fInfoLevel(-1)
  { }


  RdStationInterpolator::~RdStationInterpolator()
  { }


  fwk::VModule::ResultFlag
  RdStationInterpolator::Init()
  {
    OUT(eFinal) << "RdStationInterpolator::Init()" << endl;
    string tmpstring;
    Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdStationInterpolator");
    topBranch.GetChild("MaximumAllowedDistance").GetData(fMaximumDistance);
    topBranch.GetChild("TracePaddingAtFront").GetData(fTracePaddingAtFront);
    topBranch.GetChild("MinimumTraceLength").GetData(fMinimumTraceLength);
    topBranch.GetChild("SimEventRate").GetData(fEventRate);
    topBranch.GetChild("FirstEventTime").GetData(fFirstEventTime);
    topBranch.GetChild("UseRateTiming").GetData(tmpstring);
    fUseRateTiming = utl::StringEquivalent(tmpstring,"yes");
    topBranch.GetChild("IncludedStationIds").GetData(fIncludedStationIds);
    topBranch.GetChild("ExcludedStationIds").GetData(fExcludedStationIds);
    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);

    ostringstream info;
    if (fUseRateTiming) {
      info << "\n Ignoring event times provided by simulations! Using rate timing instead.\n"
           << " First event date " << UTCDateTime(fFirstEventTime) << "\n"
           << " Events have a rate of " << fEventRate/hertz << " Hz\n";
    }
    info << "\n Padding time series data at front by "
         << fTracePaddingAtFront/nanosecond
         << " ns and at end to a total length of at least "
         << fMinimumTraceLength/nanosecond
         << " ns."
         << endl;
    OUT(eFinal) <<info.str() << endl;

    return eSuccess;
  } //Init


  fwk::VModule::ResultFlag
  RdStationInterpolator::Run(evt::Event& theevent)
  {
    // This module should not be used to create a sim shower. The interface of the SimRadioPulse changed.
    INFO("I am a deprecated module. Please update me!");

    evt::ShowerSimData& simshow = theevent.GetSimShower();
    evt::RadioSimulation& rsim = simshow.GetRadioSimulation();

    const utl::CoordinateSystemPtr corsys = rsim.GetLocalCoordinateSystem();
    BuildSimShower(simshow,corsys);
    Detector& Det = Detector::GetInstance();
    TimeStamp eventtime = rsim.GetEventTime();

    // overwrite event time if fUseRateTiming is activated
    if (fUseRateTiming)
      eventtime = fFirstEventTime+fEventTimeOffset;

    // check if the event time has an uninitialized (illegal) value
    if (eventtime == TimeStamp(0,0)) {
      WARNING("RdStationInterpolator: Event time of RadioSimulation is not available. Aborting ...");
      return eFailure;
    }

    Det.Update(eventtime);
    const rdet::RDetector& RadioDet = Det.GetRDetector();

    if (!theevent.HasREvent())
        theevent.MakeREvent();

    revt::REvent& therev = theevent.GetREvent();

    if (!rsim.GoToFirstSimRadioPulse()) {
      WARNING("RdStationInterpolator: No SimRadioPulse in RadioSimulation. Aborting ...");
      return eFailure;
    }

    therev.GetHeader().SetId(rsim.GetEventNumber());
    therev.GetHeader().SetRunNumber(rsim.GetRunNumber());
    therev.GetHeader().SetTime(eventtime);
    simshow.MakeTimeStamp(eventtime);

    // Add Trigger Information to the Event, By default Event are assumed as externally triggered
    therev.MakeTrigger();
    therev.GetTrigger().SetExternalTrigger(true);
    // End of Trigger Information part

    // fill all Simradiopulses in one multimap
    PulseMMap mm_simpulse_x;
    PulseMMap mm_simpulse_y;
    rsim.GoToFirstSimRadioPulse();
    bool ok = true;

    do {
      const SimRadioPulse& srpulse = rsim.GetNextSimRadioPulse(ok);
      //leave the loop if there were no further SimRadioPulses
      if (!ok)
        break;

      //position of SimPulses
      double xcoord = srpulse.GetLocation().GetX(corsys);
      mm_simpulse_x.insert(PulsePair(xcoord,srpulse));
      double ycoord = (srpulse.GetLocation()).GetY(corsys);
      mm_simpulse_y.insert(PulsePair(ycoord,srpulse));
    } while (ok);

    vector<double> sPosition(3, 0.);
    vector<double> distanceBetweenSimloc(3, 0.);
    vector<double> distanceToStation(3, 0.);
    vector<revt::StationFFTDataContainer> v_simFFTDataContainer(3);
    vector<double> starttime(3, 0.);
    vector<double> binning(3, 0.);
    vector<double> areaWeight(4, 0.);

    for (rdet::RDetector::StationIterator rdIt = RadioDet.StationsBegin(); rdIt != RadioDet.StationsEnd(); ++rdIt) {
      utl::Point Position = rdIt->GetPosition();
      const int statID = rdIt->GetId();
      OUT(eObscure) << "Station " << statID;

      // if station id is listed in fExcludedStationIds, skip this simulated radio pulse (it will not be re-used for the second-closest antenna station!)
      if (!fExcludedStationIds.empty())
        if (find(fExcludedStationIds.begin(), fExcludedStationIds.end(), statID) != fExcludedStationIds.end()) {
          OUT(eObscure) << " -> not interpolated as it is in the exclude list." << endl;
          continue;
        }

      // if there is a list of fIncludedStationIds skip this simulated radio pulse if station id is not in the list
      if (!fIncludedStationIds.empty())
        if (find(fIncludedStationIds.begin(), fIncludedStationIds.end(), statID) == fIncludedStationIds.end()) {
          OUT(eObscure) << " -> station not interpolated as it is not in the include list." << endl;
          continue;
        }

      //find the nearest neighbors of simulated pulses to the station (station ID, timeseries)
      const PulseMMap& mm_pulse = FindClosestSimPulsesToStation(Position, corsys, mm_simpulse_x, mm_simpulse_y);

      if (mm_pulse.empty())
        continue;

      PulseMMap::const_iterator mmpiter = mm_pulse.begin();

      //weight of the timseries in interpolation
      areaWeight[1] = mmpiter->first;

      starttime[0] = mmpiter->second.GetStartTime();
      binning[0] = mmpiter->second.GetEfieldTimeSeries().GetBinning();
      // padding of zeros to end of timeseries
      v_simFFTDataContainer[0].GetTimeSeries() = PadTimeSeries(mmpiter->second);

      ++mmpiter;

      areaWeight[2] = mmpiter->first;
      starttime[1] = mmpiter->second.GetStartTime();
      binning[1] = mmpiter->second.GetEfieldTimeSeries().GetBinning();
      v_simFFTDataContainer[1].GetTimeSeries() = PadTimeSeries(mmpiter->second);

      ++mmpiter;

      areaWeight[3] = mmpiter->first;
      starttime[2] = mmpiter->second.GetStartTime();
      binning[2] = mmpiter->second.GetEfieldTimeSeries().GetBinning();
      v_simFFTDataContainer[2].GetTimeSeries() = PadTimeSeries(mmpiter->second);

      areaWeight[0] = areaWeight[1] + areaWeight[2] + areaWeight[3];

      // check for same binning of timeseries. Must have the same binning
      if ((fabs(binning[0] - binning[1])/binning[0]) > 1e-5 || (fabs(binning[0] - binning[2])/binning[0]) > 1e-5) {
        WARNING("RdStationInterpolator: simtimeseries don't have the same binning. Aborting ...");
        return eFailure;
      }

      //get the design frequencies from the channels and check, if for all channels the same
      double designLowerFreq = 0;
      double designLowerFreqFirstCh = 0;
      double designUpperFreq = 0;
      double designUpperFreqFirstCh = 0;
      for (rdet::Station::ChannelIterator cIt = rdIt->ChannelsBegin(); cIt!= rdIt->ChannelsEnd(); ++cIt) {
        const rdet::Channel& cDetChannel = RadioDet.GetStation(rdIt->GetId()).GetChannel(cIt->GetId());
        designLowerFreq = cDetChannel.GetDesignLowerFreq();
        designUpperFreq = cDetChannel.GetDesignUpperFreq();
        if (cIt == rdIt->ChannelsBegin()) {
          designLowerFreqFirstCh = designLowerFreq;
          designUpperFreqFirstCh = designUpperFreq;
        }
        if (designLowerFreq != designLowerFreqFirstCh || designUpperFreq != designUpperFreqFirstCh) {
          WARNING("RdStationInterpolator: design frequency different for different channels");
          return eFailure;
        }
      }

      // interpolate timeseries

      for (unsigned int i = 0; i <= v_simFFTDataContainer.size(); ++i) {
        if (v_simFFTDataContainer[i].GetNyquistZone() != 1) {
          WARNING("RdStationInterpolator: simtimeseries not in first Nyquist zone.");
          return eFailure;
        }
      }

      if (v_simFFTDataContainer[0].GetConstFrequencySpectrum().GetSize() != v_simFFTDataContainer[1].GetConstFrequencySpectrum().GetSize() ||
          v_simFFTDataContainer[0].GetConstFrequencySpectrum().GetSize() != v_simFFTDataContainer[2].GetConstFrequencySpectrum().GetSize()) {
        WARNING("RdStationInterpolator: simtimeseries don't have the same length.");
        return eFailure;
      }

      const revt::StationFFTDataContainer InterFFTDataContainer = InterpolateInTwoD(v_simFFTDataContainer, areaWeight, designLowerFreqFirstCh, designUpperFreqFirstCh);
      const revt::StationTimeSeries intertimeseries = InterFFTDataContainer.GetConstTimeSeries();

      //fill SimRadioPulse into the associated station

      therev.MakeStation(statID); // Create a new Station, already checks if the station exists
      revt::Station& Stat = therev.GetStation(statID);

      // create data structures
      Stat.MakeGPSData();     // should check before if exists?
      Stat.MakeTriggerData(); // should check before if exists?
      Stat.MakeSimData();     // should check before if exists?

      revt::StationTimeSeries& StatTimeSeries = Stat.GetStationTimeSeries() ;

      if (StatTimeSeries.GetSize()) {
        WARNING( "RDStationInterpolator: StatTimeSeries not empty!" );
        StatTimeSeries.ClearTrace();
      }

      // fill the simtimeseries into the Station TimeSeries
      StatTimeSeries.SetBinning(binning[0]);
      long numpaddedsamples = fTracePaddingAtFront/binning[0] + 0.5;

      // interpolate and set start time
      Stat.SetRawTraceStartTime(therev.GetHeader().GetTime() + TimeInterval((starttime[0]*areaWeight[1] + starttime[1]*areaWeight[2] + starttime[2]*areaWeight[3])/areaWeight[0] - numpaddedsamples*binning[0]));

      // fill with the interpolated trace
      for (revt::StationTimeSeries::ConstIterator iter = intertimeseries.Begin(); iter != intertimeseries.End(); ++iter) {
        StatTimeSeries.PushBack(*iter);
      }

    } // for (rdet::RDetector::StationIterator...)

    return eSuccess;
  } // Run


  fwk::VModule::ResultFlag
  RdStationInterpolator::Finish()
  {
    OUT(eIntermediate)<<"RdStationInterpolator::Finish()"<<endl;
    return eSuccess;
  } // Finish



  PulseMMap
  RdStationInterpolator::FindClosestSimPulsesToStation(const utl::Point& pt, const utl::CoordinateSystemPtr& coord,
                                                       const PulseMMap& mm_srpulse_x,
                                                       const PulseMMap& mm_srpulse_y)
    const
  {
    PulseMMap mm_triangle;
    PulseMMap::const_iterator iterx = mm_srpulse_x.lower_bound(pt.GetX(coord));

    //check if x value lays in simulated array
    if (iterx == mm_srpulse_x.end() || iterx == mm_srpulse_x.begin()) {
      return mm_triangle;  //return empty map
    }

    bool firsty = false;
    // find surrounding neighbors
    PulseMMap m_surroundingsrpulse = FindSurroundingNN(mm_srpulse_x, pt, coord, firsty);

    if (m_surroundingsrpulse.empty()) {
      firsty = true;
      PulseMMap::const_iterator itery = mm_srpulse_y.lower_bound(pt.GetY(coord));
      //check if y value lays in simulated array
      if (itery == mm_srpulse_y.end() || itery == mm_srpulse_y.begin()) {
        return mm_triangle; //return empty map
      }

      m_surroundingsrpulse = FindSurroundingNN(mm_srpulse_y, pt, coord, firsty);

      if (m_surroundingsrpulse.empty()) {
        return mm_triangle; //return empty map
      }
    }

    //search for surrounding triangle
    vector<utl::Point> simloc_triangle;
    vector<double> areaWeight;

    double minarea = 0;
    bool surrTrianglefound = false;
    PulseMMap::iterator it = m_surroundingsrpulse.begin();
    for (unsigned int i = 0; i < 2; ++i, ++it) {
      PulseMMap::iterator secIt = boost::next(it);
      for (unsigned int j = i + 1; j < 3; ++j, ++secIt) {
        unsigned int k = j + 1;
        for (PulseMMap::iterator thirdIt = boost::next(secIt); thirdIt != m_surroundingsrpulse.end(); ++thirdIt, ++k) {
          simloc_triangle.clear();
          simloc_triangle.push_back(it->second.GetLocation());
          simloc_triangle.push_back(secIt->second.GetLocation());
          simloc_triangle.push_back(thirdIt->second.GetLocation());
          areaWeight = GetAreaOfTriangles(pt, simloc_triangle, coord);

          if (fabs(areaWeight[1] + areaWeight[2] + areaWeight[3] - areaWeight[0]) >= 1e-5)
            continue;
          if (surrTrianglefound && minarea <= areaWeight[0])
            continue;
          minarea = areaWeight[0];

          mm_triangle.clear();
          mm_triangle.insert(PulsePair(areaWeight[2], it->second));
          mm_triangle.insert(PulsePair(areaWeight[3], secIt->second));
          mm_triangle.insert(PulsePair(areaWeight[1], thirdIt->second));
          surrTrianglefound = true;
        }
      }
    }

    return mm_triangle;
  }



  PulseMMap
  RdStationInterpolator::FindSurroundingNN(const PulseMMap& mm_simradiopulse, const utl::Point& pt,
                                           const utl::CoordinateSystemPtr& coord, const bool inverse)
    const
  {
    PulseMMap::const_iterator itera = inverse ?
      mm_simradiopulse.lower_bound(pt.GetY(coord)) :
      mm_simradiopulse.lower_bound(pt.GetX(coord));

    //both rows of simulated arrays around station a
    const double uppera = itera->first;
    --itera;
    const double lowera = itera->first;

    //get all simulated pulses with upper and lower a
    typedef PulseMMap::const_iterator PulseIterator;
    typedef pair<PulseIterator, PulseIterator> PulseIteratorPair;
    PulseIteratorPair upperb_range = mm_simradiopulse.equal_range(uppera);
    PulseIteratorPair lowerb_range = mm_simradiopulse.equal_range(lowera);

    //fill a map for each a
    PulseMMap m_upperb;
    for (PulseIterator upperb_iter = upperb_range.first;
         upperb_iter != upperb_range.second; ++upperb_iter) {
      const double b = inverse ?
        upperb_iter->second.GetLocation().GetX(coord) :
        upperb_iter->second.GetLocation().GetY(coord);
      m_upperb.insert(PulsePair(b, upperb_iter->second));
    }
    PulseMMap m_lowerb;
    for (PulseIterator lowerb_iter = lowerb_range.first;
         lowerb_iter != lowerb_range.second; ++lowerb_iter) {
      const double b = inverse ?
        lowerb_iter->second.GetLocation().GetX(coord) :
        lowerb_iter->second.GetLocation().GetY(coord);
      m_lowerb.insert(PulsePair(b, lowerb_iter->second));
    }

    //get four sourrounding nearest neighbors for station
    PulseMMap m_surroundingsrpulse;
    PulseIterator itu_b;
    PulseIterator itl_b;
    if (!inverse) {
      itu_b = m_upperb.lower_bound(pt.GetY(coord));
      itl_b = m_lowerb.lower_bound(pt.GetY(coord));
    } else {
      itu_b = m_upperb.lower_bound(pt.GetX(coord));
      itl_b = m_lowerb.lower_bound(pt.GetX(coord));
    }

    if (itu_b == m_upperb.end() || itu_b == m_upperb.begin() || itl_b == m_lowerb.end() || itl_b == m_lowerb.begin()) {
      return m_surroundingsrpulse;
    }

    m_surroundingsrpulse.insert(PulsePair(itu_b->first, itu_b->second));
    m_surroundingsrpulse.insert(PulsePair(itl_b->first, itl_b->second));
    --itu_b;
    --itl_b;
    m_surroundingsrpulse.insert(PulsePair(itu_b->first, itu_b->second));
    m_surroundingsrpulse.insert(PulsePair(itl_b->first, itl_b->second));

    return m_surroundingsrpulse;
  }



  vector<double>
  RdStationInterpolator::GetAreaOfTriangles(const utl::Point& stationpos, const vector<utl::Point>& v_simloc, const utl::CoordinateSystemPtr& coord)
    const
  {
    // distance from neighbor to station
    vector<double> distStation;
    for (int i = 0; i < 3; ++i) {
      double distx = stationpos.GetX(coord) - v_simloc[i].GetX(coord);
      double disty = stationpos.GetY(coord) - v_simloc[i].GetY(coord);
      distStation.push_back(sqrt(pow(distx, 2) + pow(disty, 2)));
    }

    // distance between neighbors
    vector<double> distSim;
    distSim.push_back((v_simloc[1] - v_simloc[0]).GetR(coord));
    distSim.push_back((v_simloc[2] - v_simloc[1]).GetR(coord));
    distSim.push_back((v_simloc[0] - v_simloc[2]).GetR(coord));

    // calculating the areas of the different triangles
    vector<double> s;
    vector<double> area;

    s.push_back((distSim[0] + distSim[1] + distSim[2])/2);
    area.push_back(sqrt(s[0] * (s[0] - distSim[0])*(s[0] - distSim[1])*(s[0] - distSim[2])));
    s.push_back((distSim[0] + distStation[0] + distStation[1])/2);
    area.push_back(sqrt(s[1] * (s[1] - distSim[0])*(s[1] - distStation[0])*(s[1] - distStation[1])));
    s.push_back((distSim[1] + distStation[1] + distStation[2])/2);
    area.push_back(sqrt(s[2] * (s[2] - distSim[1])*(s[2] - distStation[1])*(s[2] - distStation[2])));
    s.push_back((distSim[2] + distStation[2] + distStation[0])/2);
    area.push_back(sqrt(s[3] * (s[3] - distSim[2])*(s[3] - distStation[2])*(s[3] - distStation[0])));

    return area;
  }



  revt::StationTimeSeries
  RdStationInterpolator::PadTimeSeries(const SimRadioPulse& simtimeseries)
    const
  {
    revt::StationTimeSeries PaddedTimeSeries;

    Vector3D V3DZero;
    V3DZero = 0.0, 0.0, 0.0;

    double padbinning = simtimeseries.GetBinning();
    PaddedTimeSeries.SetBinning(padbinning);

    long numpaddedsamples = fTracePaddingAtFront/padbinning + 0.5;

    for (long i = 0; i < numpaddedsamples; ++i)
      PaddedTimeSeries.PushBack(V3DZero);

    for (revt::StationTimeSeries::ConstIterator iter = simtimeseries.GetEfieldTimeSeries().Begin(); iter != simtimeseries.GetEfieldTimeSeries().End(); ++iter) {
      PaddedTimeSeries.PushBack(*iter);
    }

    // now pad the station trace at the end up to the given minimum trace length
    while (fMinimumTraceLength - PaddedTimeSeries.GetSize()*padbinning > 1e-6)
      PaddedTimeSeries.PushBack(V3DZero);

    // if number of samples is uneven after padding, pad one more sample (otherwise it will again be clipped at the next FFT)
    if (PaddedTimeSeries.GetSize()/2*2 != PaddedTimeSeries.GetSize())
      PaddedTimeSeries.PushBack(V3DZero);

    return PaddedTimeSeries;
  }



  revt::StationFFTDataContainer
  RdStationInterpolator::InterpolateInTwoD(const std::vector<revt::StationFFTDataContainer>& simData,
                                           const std::vector<double>& area, const double designlowerfreq, const double designupperfreq)
    const
  {
    const double specbinningOne = simData[0].GetConstFrequencySpectrum().GetBinning();
    const StationFrequencySpectrum::SizeType size1 = simData[0].GetConstFrequencySpectrum().GetSize();
    const StationFrequencySpectrum& freq1 = simData[0].GetConstFrequencySpectrum();

    // unused    const StationFrequencySpectrum::SizeType size2 = simData[1].GetConstFrequencySpectrum().GetSize();
    const StationFrequencySpectrum& freq2 = simData[1].GetConstFrequencySpectrum();

    // unused    const StationFrequencySpectrum::SizeType size3 = simData[2].GetConstFrequencySpectrum().GetSize();
    const StationFrequencySpectrum& freq3 = simData[2].GetConstFrequencySpectrum();

    TraceV3C V_interspectrum;
    V_interspectrum.SetBinning(specbinningOne);
    Vector3C element;

    Vector3D phaseElementOne;
    Vector3D phaseElementTwo;
    Vector3D phaseElementThree;
    Vector3D phaseElementInter;

    Vector3D jumps1;
    jumps1 = 0.0, 0.0, 0.0;
    Vector3D jumps2;
    jumps2 = 0.0, 0.0, 0.0;
    Vector3D jumps3;
    jumps3 = 0.0, 0.0, 0.0;

    for (StationFrequencySpectrum::SizeType i = 0 ; i < size1 ; ++i) {
      for(unsigned int j = 0; j < 3; ++j) { // loop over x,y and z value
        if (i >= designlowerfreq/specbinningOne && i <= designupperfreq/specbinningOne) {
          // interpolate amplitudes
          double interamplitude = (abs(freq1[i][j])*area[1] + abs(freq2[i][j])*area[2] + abs(freq3[i][j])*area[3])/area[0];

          // unwrap phases of initial signals
          if (i > 1) {

            if (arg(freq1[i][j]) < arg(freq1[i-1][j]) && jumps1[j] >= 0) {
              if (fabs(arg(freq1[i][j])-arg(freq1[i-1][j])) > fabs(arg(freq1[i][j]) + utl::kPi - arg(freq1[i-1][j]))) {
                ++jumps1[j];
              }
            } else {
              if (fabs(arg(freq1[i][j])-arg(freq1[i-1][j])) > fabs(arg(freq1[i][j]) - utl::kPi - arg(freq1[i-1][j]))) {
                --jumps1[j];
              }
            }

            if (arg(freq2[i][j]) < arg(freq2[i-1][j]) && jumps2[j] >= 0) {
              if (fabs(arg(freq2[i][j]) - arg(freq2[i-1][j])) > fabs(arg(freq2[i][j]) + utl::kPi - arg(freq2[i-1][j]))) {
                ++jumps2[j];
              }
            } else {
              if (fabs(arg(freq2[i][j]) - arg(freq2[i-1][j])) > fabs(arg(freq2[i][j]) - utl::kPi - arg(freq2[i-1][j]))) {
                --jumps2[j];
              }
            }

            if (arg(freq3[i][j]) < arg(freq3[i-1][j]) && jumps3[j] >= 0) {
              if (fabs(arg(freq3[i][j]) - arg(freq3[i-1][j])) > fabs(arg(freq3[i][j]) + utl::kPi - arg(freq3[i-1][j]))) {
                ++jumps3[j];
              }
            } else {
              if (fabs(arg(freq3[i][j]) - arg(freq3[i-1][j])) > fabs(arg(freq3[i][j]) - utl::kPi - arg(freq3[i-1][j]))) {
                --jumps3[j];
              }
            }

          }

          phaseElementOne[j] = arg(freq1[i][j]) + jumps1[j]*2*utl::kPi;
          phaseElementTwo[j] = arg(freq2[i][j]) + jumps2[j]*2*utl::kPi;
          phaseElementThree[j] = arg(freq3[i][j]) + jumps3[j]*2*utl::kPi;

          // interpolate phases
          phaseElementInter[j] = (phaseElementOne[j]*area[1]  + phaseElementTwo[j]*area[2] + phaseElementThree[j]*area[3])/area[0];

          // interpolated spectrum
          element[j] = polar(interamplitude, phaseElementInter[j]);

        } else {
          element[j] = 0;
        }
      }

      V_interspectrum.PushBack(element);
    }

    revt::StationFFTDataContainer interpolatedFFTDataContainer;
    interpolatedFFTDataContainer.GetFrequencySpectrum() = V_interspectrum;

    return interpolatedFFTDataContainer;
  }



  void
  RdStationInterpolator::BuildSimShower(evt::ShowerSimData& theshower, const utl::CoordinateSystemPtr& coord)
  {
    // This function will build a Sim shower with correct data
    // Inspired by the EventGenerator Module

    if (!theshower.HasPosition())
      theshower.MakeGeometry(Point(0, 0, 0, coord));

  }

} // nameSpace