RdChannelBeaconTimingCalibrator.cc

Go to the documentation of this file.

#include "RdChannelBeaconTimingCalibrator.h"

#include <fwk/CentralConfig.h>

#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/config.h>
#include <utl/AugerUnits.h>
#include <utl/MathConstants.h>
#include <utl/TraceAlgorithm.h>
#include <utl/TimeStamp.h>
#include <utl/TimeInterval.h>
#include <utl/FFTDataContainerAlgorithm.h>
#include <utl/CoordinateSystem.h>
#include <utl/PhysicalConstants.h>
#include <utl/Point.h>
#include <utl/Vector.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 <det/Detector.h>
#include <rdet/RDetector.h>

#include <rdet/Station.h>
#include <rdet/Channel.h>

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

#include <CLHEP/Matrix/Matrix.h>
#include <CLHEP/Matrix/Vector.h>

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

const double speedOfLight = kSpeedOfLight / 1.00025; // for adjusting refractive index

namespace RdChannelBeaconTimingCalibrator {

  VModule::ResultFlag
  RdChannelBeaconTimingCalibrator::Init()
  {
    const Branch topBranch =
      CentralConfig::GetInstance()->GetTopBranch("RdChannelBeaconTimingCalibrator");

    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
    topBranch.GetChild("BeaconFrequencies").GetData(fBeaconFrequencies);
    topBranch.GetChild("ReferenceChannel").GetData(fReferenceChannel);
    topBranch.GetChild("FixedReferenceStationID").GetData(fFixedReferenceStationID);
    topBranch.GetChild("NoiseFilterSize").GetData(fNoiseFilterSize);
    topBranch.GetChild("MinSNR").GetData(fMinSNR);
    topBranch.GetChild("MaxTimeDifference").GetData(fMaxTimeDifference);
    topBranch.GetChild("Epsilon").GetData(fEpsilon);
    topBranch.GetChild("MinGoodFreqs").GetData(fMinGoodFreqs);
    topBranch.GetChild("EqualizeStationTimeStamps").GetData(fEqualizeStationTimeStamps);
    topBranch.GetChild("InputPhaseFile").GetData(fInputPhaseFile);
    // read in the debug station list as vector, and then copy it to the set class variable
    vector<int> dummyStationIDlist;
    topBranch.GetChild("DebugStationIds").GetData(dummyStationIDlist);
    fDebugStationIds = set<int>(dummyStationIDlist.begin(), dummyStationIDlist.end());
    topBranch.GetChild("DebugOutputFile").GetData(fDebugOutputFile);
    topBranch.GetChild("OverwriteDebugFile").GetData(fOverwriteDebugFile);
    topBranch.GetChild("RejectIfNoRefPhaseAvailable").GetData(fRejectIfNoRefPhaseAvailable);
    topBranch.GetChild("RejectIfCorrectionFailed").GetData(fRejectIfCorrectionFailed);
    topBranch.GetChild("UseCrossCorrelationMethod").GetData(fUseCrossCorrelationMethod);
    topBranch.GetChild("UseBeaconDistanceInsteadOfDatabaseValues").GetData(fUseBeaconDistanceInsteadOfDatabaseValues);
    topBranch.GetChild("BeaconPosition").GetData(fBeaconPosition);
    topBranch.GetChild("BeaconTimeCalibrationUncertainty").GetData(fBeaconTimeCalibrationUncertainty);
    topBranch.GetChild("StartingDateTimeForBeaconCorrection").GetData(fStartingDateTimeForBeaconCorrection);

    //check if the  beacon position has the right format
    if(fBeaconPosition.size() < 3) {
      ERROR("Incorrect definition of beacon position in XML file: One or more coordinates are missing!");
      return eBreakLoop;
    }

    return eSuccess;
  }


  VModule::ResultFlag
  RdChannelBeaconTimingCalibrator::Run(evt::Event& event)
  {
    // variables needed during the whole method

    std::map<int, std::map<int, double>> phases;

    std::map<int, std::map<int, double>> modPhases;

    std::map<int, std::map<int, double>> timeDiffsFreq;

    std::map<int, std::map<int, double>> SNRs;

    std::map<int, double> timeOffsets;

    std::map<int, double> timeOffsetUncertainty;

    // Check if there are events at all
    if (!event.HasREvent()) {
      WARNING("RdChannelBeaconTimingCalibrator::No radio event found!");
      return eContinueLoop;
    } else {
      REvent& rEvent = event.GetREvent();
      stringstream message;
      message << "Radio event found with "
               << rEvent.GetNumberOfStations()
               << " stations!";
      if (fInfoLevel >= eInfoDebug)
        INFO(message.str());

      // Read event header and check event time
      const Header& rHeader = rEvent.GetHeader();
      message.str("");
      message << "Found header with ID "
               << rHeader.GetId()
               << " and timestamp: "
               << rHeader.GetTime();
      if (fInfoLevel >= eInfoDebug)
        INFO(message.str());

      if (rHeader.GetTime() < fStartingDateTimeForBeaconCorrection) {
        message.str("");
        message << "Skipped event with ID "
                << rHeader.GetId()
                << " because event time is before: "
                << fStartingDateTimeForBeaconCorrection;
        if (fInfoLevel >= eInfoFinal)
          INFO(message.str());
        return eSuccess;
      }

      // add statistics of beacon successes per stationID
      // at the beginning assume a failure and then substract it at the end, if correction was successful
      for (REvent::ConstStationIterator sIt = rEvent.StationsBegin();
           sIt != rEvent.StationsEnd(); ++sIt) {
        const Station& station = *sIt;
        // check, if station is already tracked, and if not intialize its statistics counter
        if (fBeaconCorrectionSuccesses.find(station.GetId()) == fBeaconCorrectionSuccesses.end() )
          fBeaconCorrectionSuccesses[station.GetId()] = 0;
        if (fBeaconCorrectionFailures.find(station.GetId()) == fBeaconCorrectionFailures.end() )
          fBeaconCorrectionFailures[station.GetId()] = 0;
        // assume failure (and substract it later, if correction successful)
        ++fBeaconCorrectionFailures[station.GetId()];
      }

      // get Detector
      const det::Detector& Det = det::Detector::GetInstance();
      const rdet::RDetector& RDetector = Det.GetRDetector();

      // 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
      std::vector<double> beaconFrequencies(fBeaconFrequencies);
      if(beaconFrequencies.size() < 1) {
        try {
          beaconFrequencies = RDetector.GetBeaconFrequencies(); // for example for 2013 returns: 0.058887 0.071191 0.068555 0.061523
        } catch (utl::DataNotFoundInDBException& err) {
          ERROR ("Error reading Beacon frequencies from database, skipping event!");
          return eContinueLoop;
        }
      }

      // read the reference phases (from txt file or database)
      std::vector<std::map<int, double> > referencePhases;
      if (fInputPhaseFile != "") {
        referencePhases = getReferencePhases(beaconFrequencies);
      } else {
        referencePhases = vector<map<int, double> > (beaconFrequencies.size());
        for (REvent::ConstStationIterator sIt = rEvent.StationsBegin();
          sIt != rEvent.StationsEnd(); ++sIt) {
          const Station& station = *sIt;
          const int stationId = station.GetId();
          try {
            for (unsigned int beaconFreqN=0; beaconFreqN < beaconFrequencies.size(); ++beaconFreqN) {
              referencePhases[beaconFreqN][stationId] = RDetector.GetBeaconReferencePhase(stationId, beaconFrequencies[beaconFreqN]);
            }
          } catch (utl::DataNotFoundInDBException& err) {
            message.str("");
            message << "Could not read beacon reference phases for station "
                    << station.GetId()
                    << " from database.";
           if (fInfoLevel >= eInfoFinal)
             WARNING(message.str());
          }
        }
      }
      // final cross check
      if (referencePhases.size() != beaconFrequencies.size()) {
        ERROR("Could not read reference phases correctly.");
        return eFailure;
      }

      // check if there are reference phases for at least two stations
      if (referencePhases[0].size() < 2) {
        message.str("");
        message << "Event has only "
                << referencePhases[0].size()
                << " stations with available reference phases. At least 2 are needed for the beacon correction.";
        ERROR(message.str());
        return eContinueLoop;
      }

      // check if fixed reference station is contained in the event
      if (fFixedReferenceStationID > 0)
        if (!rEvent.HasStation(fFixedReferenceStationID)) {
          message.str("");
          message << "Reference station with ID "
               << fFixedReferenceStationID
               << " is not contained in the event: No Beacon correction possible!";
          ERROR(message.str());
          return eContinueLoop;
        }

      // vector for beacon frequencies acutally corresponding to the
      vector<double> beaconFreqsInSpectrum;  // attention: can change from event to event, since trace length is variable

      // loop through stations and for every station through every channel
      // calculate the phase differences for the reference channel (default: NS, high gain = channel 2)
      unsigned int numberOfNotRejectedStations = 0;
      for (auto sIt = rEvent.StationsBegin(); sIt != rEvent.StationsEnd(); ++sIt) {
        Station& station = *sIt;

        // check if all references phases are available for this station
        bool skipStation = false;
        for (unsigned int beaconFreqN = 0; beaconFreqN < beaconFrequencies.size(); ++beaconFreqN) {
          if (referencePhases[beaconFreqN].find(station.GetId()) == referencePhases[beaconFreqN].end() ) {
            message.str("");
            message << "No reference phase for station "
                    << station.GetId()
                    << " and frequency #"
                    << beaconFreqN+1
                    << " available!";
                    if (fRejectIfNoRefPhaseAvailable)
                      message << " Station will be rejected.\n";
                    else
                      message << ".\n";
            if (fInfoLevel >= eInfoFinal)
              WARNING(message.str());
            if (fRejectIfNoRefPhaseAvailable)
              station.SetRejectedReason(eBeaconCorrectionFailed);
            skipStation = true;
            break;
          }

          // check if beacon frequencies are inside of frequency spectrum of the station
          // (which is not the case for low frequency stations)
          if (    (beaconFrequencies[beaconFreqN] < RDetector.GetStation(station.GetId()).GetChannel(fReferenceChannel).GetDesignLowerFreq())
               || (beaconFrequencies[beaconFreqN] > RDetector.GetStation(station.GetId()).GetChannel(fReferenceChannel).GetDesignUpperFreq()) ) {
            message.str("");
            message << "Beacon frequency "
                    << beaconFrequencies[beaconFreqN]/megahertz
                    << " MHz is outside of the design frequency band of station "
                    << station.GetId()
                    << ". --> Station will be rejected.\n";
            if (fInfoLevel >= eInfoFinal)
              WARNING(message.str());
            station.SetRejectedReason(eBeaconCorrectionFailed);
            skipStation = true;
            break;
          }
        }
        if (skipStation)
          continue;

        // Get spectrum of reference channel
        Channel& channel = station.GetChannel(fReferenceChannel);

        message.str("");
        message << "Calculating the phases of the beacon frequencies for channel "
                << channel.GetId()
                << " of station "
                << channel.GetStationId()
                << ".";
        if (fInfoLevel >= eInfoDebug)
          INFO(message.str());

        const auto spectrum = channel.GetChannelFrequencySpectrum();

        // check if spectrum has the correct length, such that the frequency bins are at the expected values
        // for which the beacon has been optimized
        // are more general implementation might follow in the future
        if ( ((spectrum.GetSize()-1) % 1024) != 0 ) {
          message.str("");
          message << "Trace should consist of 2048 samples or an integer multiple of it! It has:"
                  << (spectrum.GetSize()-1)*2;
          if (fInfoLevel >= eInfoFinal)
            WARNING(message.str());
          //return eFailure; // lets try even if frequency bins are non-optimal
          continue; // just ignore stations with improper trace lengths (e.g., Dutch stations with 1024 samples, only)
        }

        // get the amplitude spectrum for SNR calculation
        Trace<double> amplitudeSpectrum;
        for (ChannelFrequencySpectrum::SizeType i = 0; i < spectrum.GetSize(); ++i)
          amplitudeSpectrum.PushBack(abs(spectrum[i]));


        // loop through beacon frequencies, find corresponding frequency in the spectrum
        // and store phases
        for (unsigned int beaconFreqN=0; beaconFreqN < beaconFrequencies.size(); ++beaconFreqN) {
          // check if frequency is inside the range (take care that the order of the spectrum is unknown)
          if(  (beaconFrequencies[beaconFreqN] <
                min(channel.GetFrequencyOfBin(0),channel.GetFrequencyOfBin(spectrum.GetSize()-1)))
            || (beaconFrequencies[beaconFreqN] >
                max(channel.GetFrequencyOfBin(0),channel.GetFrequencyOfBin(spectrum.GetSize()-1))) ) {
            ERROR("Beacon frequency must be inside of the spectrum!");
            return eFailure;
          }
          // find nearest frequency bin:
          double bin = (beaconFrequencies[beaconFreqN]-channel.GetFrequencyOfBin(0)) * (spectrum.GetSize()-1)
                    / (channel.GetFrequencyOfBin(spectrum.GetSize()-1) - channel.GetFrequencyOfBin(0));

          // store frequencies of the spectrum bins closest to the real beacon frequencies
          // (needed for later conversion of phase differences to time differences)
          // only store frequencies for the first station (assuming that all stations have the same trace lengths)
          if (beaconFreqsInSpectrum.size() < beaconFrequencies.size())
            beaconFreqsInSpectrum.push_back( channel.GetFrequencyOfBin(round(bin)) );

          // consistency check whether frequency of the spectrum bin is still the same (within 1e-6 due to double imprecision)
          if ( abs(beaconFreqsInSpectrum[beaconFreqN] / channel.GetFrequencyOfBin(round(bin)) - 1.) > 1e-6) {
            message.str("");
            message << "Frequency of bin corresponding to beacon frequency "
                    << beaconFreqN + 1
                    << " changed from "
                    << beaconFreqsInSpectrum[beaconFreqN]/megahertz << "MHz"
                    << " to "
                    << channel.GetFrequencyOfBin(round(bin))/megahertz << "MHz";
            ERROR(message.str());
            return eFailure;
          }


          // get the phase, amplitude and noise level of the bin
          double phase = arg(spectrum[round(bin)]);
          double amplitude = abs(spectrum[round(bin)]);

          // substract reference phase
          double modPhase = phase - referencePhases[beaconFreqN][channel.GetStationId()];
          if (modPhase > kPi) {
            modPhase -= kTwoPi;
          } else if (modPhase < -kPi) {
            modPhase += kTwoPi;
          }

          //Store original and modified phases
          phases[beaconFreqN][station.GetId()] = phase;
          modPhases[beaconFreqN][station.GetId()] = modPhase;

          //Calculate and store SNR
          double noise = TraceAlgorithm::Median(amplitudeSpectrum,
                        round(bin)-trunc(fNoiseFilterSize/2.0),
                        round(bin)+ceil(fNoiseFilterSize/2.0)-1.);
          double snr = (amplitude*amplitude)/(noise*noise) ;
          SNRs[beaconFreqN][station.GetId()] = snr;


          // store SNRs for debugging, monitoring etc.
          if (fBeaconSNRs.find(station.GetId()) == fBeaconSNRs.end())
            fBeaconSNRs[station.GetId()] = std::map<int,double>();

          if (fBeaconSNRs[station.GetId()].find(beaconFreqN) == fBeaconSNRs[station.GetId()].end())
            fBeaconSNRs[station.GetId()][beaconFreqN] = snr;
          else
            fBeaconSNRs[station.GetId()][beaconFreqN] += snr;

        }
        ++numberOfNotRejectedStations;

        // increase counter for SNR normalization
        if (fBeaconSNRcounter.find(station.GetId()) == fBeaconSNRcounter.end() )
          fBeaconSNRcounter[station.GetId()] = 1;
        else
          ++fBeaconSNRcounter[station.GetId()];
      }

      // check if there are at least two not rejected stations
      if (numberOfNotRejectedStations < 2) {
        message.str("");
        message << "Event has only "
                << numberOfNotRejectedStations
                << " stations which are not rejected. At least 2 are needed for the beacon correction.";
        ERROR(message.str());
        return eContinueLoop;
      }

      // find reference station and get reference phases:
      // Either take largest SNR in beacon frequency with lowest SNR, or use pre-defined fixed reference station
      double bestSNR = 0;
      int refStationID = 0;
      double refStationTimeStamp = 0;    // GPS time (relative to event time) of reference station
      double refStationPhases[beaconFrequencies.size()];

      for (REvent::StationIterator sIt = rEvent.StationsBegin();
           sIt != rEvent.StationsEnd(); ++sIt) {
        Station& station = *sIt;

        if ((fFixedReferenceStationID>0) && (station.GetId()!=fFixedReferenceStationID))
          continue;

        // skip stations without beacon reference phases
        bool skipStation = false;
        for (unsigned int beaconFreqN=0; beaconFreqN < beaconFrequencies.size(); ++beaconFreqN) {
          if ( referencePhases[beaconFreqN].find(station.GetId()) == referencePhases[beaconFreqN].end() )
            skipStation = true;
        }
        if (skipStation) {
          //cout << "Station " << station.GetId() << " skipped because of missing reference phases!" << endl;
          continue;
        }

        // find beacon frequency with smallest SNR
        double smallestSNRofFreq = 1e99;

        for (unsigned int beaconFreqN=0; beaconFreqN < beaconFrequencies.size(); ++beaconFreqN) {
          if (SNRs[beaconFreqN][station.GetId()] < smallestSNRofFreq)
            smallestSNRofFreq = SNRs[beaconFreqN][station.GetId()];
        }

        // if this is the best station so far, than store reference information
        if (smallestSNRofFreq > bestSNR) {
          bestSNR = smallestSNRofFreq;
          refStationID = station.GetId();
          refStationTimeStamp = station.GetRecData().GetParameter(eTraceStartTime);
          for (unsigned int beaconFreqN=0; beaconFreqN < beaconFrequencies.size(); ++beaconFreqN)
            refStationPhases[beaconFreqN] = modPhases[beaconFreqN][station.GetId()];
        }


        // Message for debugging
        message.str("");
        message << "SNR (in power) for station "
                << station.GetId()
                << " : ";
        for (unsigned int beaconFreqN=0; beaconFreqN < beaconFrequencies.size(); ++beaconFreqN)
          message << SNRs[beaconFreqN][station.GetId()] << " \t ";
        if (fInfoLevel >= eInfoDebug)
          INFO(message.str());
      }

      // Message for debugging
      message.str("");
      message << "Selected station "
              << refStationID
              << " as reference, since its SNR (power) in the worst frequency is: "
              << bestSNR;
      if ((fInfoLevel >= eInfoIntermediate) && (fFixedReferenceStationID>0))
        INFO(message.str());

      // check if SNR is sufficient
      if (bestSNR < fMinSNR) {
        message.str("");
        message << "The SNR (power) of the worst frequency of the reference station is insufficient: "
                << bestSNR
                << " (should be >= "
                << fMinSNR
                << ")!";
        if (fInfoLevel >= eInfoFinal)
          WARNING(message.str());
      }


      // get positions and expectations for delay
      double referenceDelay = 0; // only important if expectation from distances is used instead of reference phases

      if (fUseBeaconDistanceInsteadOfDatabaseValues) {

        const Point refStationPosition(RDetector.GetStation(refStationID).GetPosition());

        const Point beaconPosition(fBeaconPosition[0], fBeaconPosition[1], fBeaconPosition[2], Det.GetReferenceCoordinateSystem());

        //calculation of the referenceDelay
        referenceDelay = (refStationPosition - beaconPosition).GetMag() / speedOfLight;

        //cout << "\nTime to reference station: " << referenceDelay << " ns.\n" << endl;
      }

      // calculate phase differences and time differences
      for (REvent::StationIterator sIt = rEvent.StationsBegin();
           sIt != rEvent.StationsEnd(); ++sIt) {
        Station& station = *sIt;

        timeOffsets[station.GetId()] = 0;   // Reset time difference of this station
        timeOffsetUncertainty[station.GetId()] = 99.;   // Reset uncertainty

        // skip stations without beacon reference phases
        bool skipStation = false;
        for (unsigned int beaconFreqN=0; beaconFreqN < beaconFrequencies.size(); ++beaconFreqN) {
          if ( referencePhases[beaconFreqN].find(station.GetId()) == referencePhases[beaconFreqN].end() )
            skipStation = true;
        }
        if (skipStation) {
          //cout << "Station " << station.GetId() << " skipped because of missing reference phases!" << endl;
          continue;
        }

        // only calculate time difference if at least 2 beacon frequencies with good SNR are available
        // also check if reference phases are available for station
        unsigned int goodSNRfreqs = 0;

        for (unsigned int beaconFreqN=0; beaconFreqN < beaconFrequencies.size(); ++beaconFreqN) {
          if ((SNRs[beaconFreqN][station.GetId()] > fMinSNR) && (SNRs[beaconFreqN][refStationID] > fMinSNR))
            ++goodSNRfreqs;
        }

        if (goodSNRfreqs < fMinGoodFreqs) {
          message.str("");
          message << "Station " << station.GetId() << " has " << goodSNRfreqs
                  << " beacon signals with good SNR. At least "
                  << fMinGoodFreqs << " are requiered."
                  << " No timing correction possible.";
          if (fInfoLevel >= eInfoFinal)
            WARNING(message.str());
          // set time offset for this station to 0, i.e. perform no timing correction
          timeOffsets[station.GetId()] = 0;
          continue; // next station
        }

        // Take possible difference of GPS time stamps into account
        double timeStampOffsetGPS = station.GetRecData().GetParameter(eTraceStartTime) - refStationTimeStamp;  // ns = Auger base unit

        // in rare cases the GPS time stamp of the station does not match the GPS time stamp of the event
        // --> ignore those stations, if mismatch is larger than 0.01 seconds (corresponds to 3000 km light travel)
        if (abs(timeStampOffsetGPS) > 1e7) {
          message.str("");
          message << "Trace start time of station "
                  << station.GetId()
                  << " is "
                  << timeStampOffsetGPS
                  << " ns larger than event time stamp! --> Station ignored!";
          if (fInfoLevel >= eInfoFinal)
            WARNING(message.str());
          continue;
        }


        // semi-analytically find maximum of cross-correlation
        // this is the shift tau, for which the following sum is maximum:
        // sum of [freq_i * cos(phasediff_i - omega_i*tau)]

        if (fUseCrossCorrelationMethod) {
          double expectedDelay = 0; // only important if expectation from distances is used instead of reference phases

          // get delay to station and relative to reference stationID
          if (fUseBeaconDistanceInsteadOfDatabaseValues) {
            const Point stationPosition(RDetector.GetStation(station.GetId()).GetPosition());

            const Point beaconPosition(fBeaconPosition[0], fBeaconPosition[1], fBeaconPosition[2], Det.GetReferenceCoordinateSystem());

            expectedDelay = (stationPosition - beaconPosition).GetMag() / speedOfLight - referenceDelay;
          }
          //cout << "\nExpected delay for station " << station.GetId() << " (relative to RefStation): " << expectedDelay << " ns." << endl;

          // search for maximum in different step sizes
          const double mainstepsize = 0.5; // stepsize in ns for rough search of maxima
          const double finestepsize = 0.01; // stepsize in ns for fine determination of maxima
          double maxShiftValue = 0;
          double lastSum = -1e99;
          double secondLastSum = -1e99;
          double currentMaxSum = -1e99;
          for (double shift = - fMaxTimeDifference; shift <= fMaxTimeDifference; shift += mainstepsize) {
            double currentSum = 0;
            for (unsigned int beaconFreqN=0; beaconFreqN < beaconFrequencies.size(); ++beaconFreqN) {
              double phasediff = phases[beaconFreqN][station.GetId()] - phases[beaconFreqN][refStationID]
                                 - timeStampOffsetGPS * (kTwoPi * beaconFreqsInSpectrum[beaconFreqN])
                                 - shift * (kTwoPi * beaconFreqsInSpectrum[beaconFreqN]);
              if (fUseBeaconDistanceInsteadOfDatabaseValues)
                phasediff += expectedDelay * (kTwoPi * beaconFreqsInSpectrum[beaconFreqN]);
              else
                phasediff += referencePhases[beaconFreqN][refStationID] - referencePhases[beaconFreqN][station.GetId()];

              // use phases calculated by distance to beacon
              // double phasediff = phases[beaconFreqN][station.GetId()] - phases[beaconFreqN][refStationID]
              //                  - timeStampOffsetGPS * (kTwoPi * beaconFreqsInSpectrum[beaconFreqN])
              //                  - (shift-expectedDelay) * (kTwoPi * beaconFreqsInSpectrum[beaconFreqN]);
              // use reference phase from data base
              // double phasediff =   phases[beaconFreqN][station.GetId()] - referencePhases[beaconFreqN][station.GetId()]
              //                   - phases[beaconFreqN][refStationID] + referencePhases[beaconFreqN][refStationID]
              //                   - timeStampOffsetGPS * (kTwoPi * beaconFreqsInSpectrum[beaconFreqN])
              //                   - shift * (kTwoPi * beaconFreqsInSpectrum[beaconFreqN]);

              // frequency factor comes from analytic calculation of cross-correlation integral:
              // --> try without this factor, to give all frequencies equal weight
              // instead reduce weight for frequencies with small SNR (for station and reference station
              double weight = 1.;
              if (SNRs[beaconFreqN][station.GetId()] < fMinSNR)
                weight *= sqrt(SNRs[beaconFreqN][station.GetId()]/fMinSNR);
              if (SNRs[beaconFreqN][refStationID] < fMinSNR)
                weight *= sqrt(SNRs[beaconFreqN][refStationID]/fMinSNR);
              currentSum += weight*cos(phasediff);// * beaconFreqsInSpectrum[beaconFreqN];

              //cout << "Sum for shift of " << shift << " ns: " << currentSum << endl;
            }

            // look if local maximum was found, then do fine search
            if ((lastSum>currentSum)&&(lastSum>secondLastSum)) {
              for (double fineshift = shift-2*mainstepsize; fineshift < shift; fineshift += finestepsize) {
                double newFineSearchSum = 0;
                for (unsigned int beaconFreqN=0; beaconFreqN < beaconFrequencies.size(); ++beaconFreqN) {
                  double phasediff = phases[beaconFreqN][station.GetId()] - phases[beaconFreqN][refStationID]
                                    - timeStampOffsetGPS * (kTwoPi * beaconFreqsInSpectrum[beaconFreqN])
                                    - fineshift * (kTwoPi * beaconFreqsInSpectrum[beaconFreqN]);
                  if (fUseBeaconDistanceInsteadOfDatabaseValues)
                    phasediff += expectedDelay * (kTwoPi * beaconFreqsInSpectrum[beaconFreqN]);
                  else
                    phasediff += referencePhases[beaconFreqN][refStationID] - referencePhases[beaconFreqN][station.GetId()];

                  double weight = 1.;
                  if (SNRs[beaconFreqN][station.GetId()] < fMinSNR)
                    weight *= sqrt(SNRs[beaconFreqN][station.GetId()]/fMinSNR);
                  if (SNRs[beaconFreqN][refStationID] < fMinSNR)
                    weight *= sqrt(SNRs[beaconFreqN][refStationID]/fMinSNR);
                  newFineSearchSum += weight*cos(phasediff);
                }

                if (newFineSearchSum > currentMaxSum) {
                  maxShiftValue = fineshift;
                  currentMaxSum = newFineSearchSum;
                }
              }
            } // if ((lastSum>currentSum)&&(lastSum>secondLastSum))

            // save old values for last to steps
            secondLastSum = lastSum;
            lastSum = currentSum;
          } // for (double shift = - fMaxTimeDifference; shift <= fMaxTimeDifference; shift += mainstepsize)
          timeOffsets[station.GetId()] = maxShiftValue;
          timeOffsetUncertainty[station.GetId()] = 0;

        } else {       //  if (fUseCrossCorrelationMethod)

          // Determine time difference, starting with a time offset of 0
          double timeOffset = 0;

          // search as long as the time offset is in the allowed range
          // while loop is quit by "break" if correct time difference is found
          while ( abs(timeOffset) <= fMaxTimeDifference) {

            // calculate time difference for each beacon frequency
            for (unsigned int beaconFreqN=0; beaconFreqN < beaconFrequencies.size(); ++beaconFreqN) {
              double phasediff = modPhases[beaconFreqN][station.GetId()] - refStationPhases[beaconFreqN]
                                - timeOffset * (kTwoPi * beaconFreqsInSpectrum[beaconFreqN])
                                - timeStampOffsetGPS * (kTwoPi * beaconFreqsInSpectrum[beaconFreqN]);

              // due to GPS time stamp offsets, phase differences can be much larger than +/- pi
              while (phasediff > kPi)
                phasediff -= kTwoPi;
              while (phasediff < -kPi)
                phasediff += kTwoPi;

              double timediff = phasediff / ( kTwoPi * beaconFreqsInSpectrum[beaconFreqN] ) + timeOffset;
              timeDiffsFreq[beaconFreqN][station.GetId()] = timediff;
            }

            // put all usable time differences for this station in one vector
            vector<double> stationTimeDiffs;
            for (unsigned int beaconFreqN=0; beaconFreqN < beaconFrequencies.size(); ++beaconFreqN)
              if ((SNRs[beaconFreqN][station.GetId()] > fMinSNR) && (SNRs[beaconFreqN][refStationID] > fMinSNR))
                stationTimeDiffs.push_back(timeDiffsFreq[beaconFreqN][station.GetId()]);


            // Debug output write time differences for this station
            message.str("");
            message << "Timediffs for station "
                    << station.GetId()
                    << " [ns]: ";
            for (unsigned int beaconFreqN=0; beaconFreqN < beaconFrequencies.size(); ++beaconFreqN)
              message << timeDiffsFreq[beaconFreqN][station.GetId()] << " \t";
            if (fInfoLevel >= eInfoDebug)
              INFO(message.str());

            // test if all time differences with valid SNR are within the acceptable range (epsilon)
            // and calculate absolute mean deviation
            double meanTdiff = accumulate(stationTimeDiffs.begin(), stationTimeDiffs.end(), 0.) / double(stationTimeDiffs.size());
            double deviationFromMean = 0;
            // message for debugging
            message.str("");
            message << "Mean tdiff = "
                    << meanTdiff
                    << " ns, deviations per beacon freq: ";

            bool tDiffsGood = true;
            for (vector<double>::iterator it = stationTimeDiffs.begin();
                      it != stationTimeDiffs.end(); ++it) {
              message << abs(*it - meanTdiff)
                      << " \t";
              deviationFromMean += abs(*it - meanTdiff) / stationTimeDiffs.size();

              // test if all time differences agree within the allowed range (epsilon)
              if (abs(*it - meanTdiff) > fEpsilon)
                tDiffsGood = false;
            }
            if (fInfoLevel >= eInfoDebug) {
              INFO(message.str());
              message.str("");
              message << "Absolute mean deviation = "
                      << deviationFromMean
                      << " ns.";
              INFO(message.str());
            }

            // If all time differences agree within the allowed range (epsilon), we are done
            // if not, try a different time offset (alternating + and -), and try again, if
            // the resulting time differences at each beacon frequency are consistent within epsilon.
            if(tDiffsGood) {
              timeOffsets[station.GetId()] = meanTdiff;
              timeOffsetUncertainty[station.GetId()] = deviationFromMean;
              break;
            } else {
              if (timeOffset > 0.) {
                timeOffset *= -1.;
              } else {
                timeOffset *= -1.;
                timeOffset += 12.*nanosecond;
                // TODO: is 12 ns a useful step size, can it be larger?
                // assumption: as long as it is shorter than the period of the highest frequency, it should be ok
                // 80 MHz --> 12.5 ns, thus 12 ns should be ok, but probably it could be even larger
              }
            }
          } // end of while
        } // else (fUseCrossCorrelationMethod)
      } // for (REvent::StationIterator sIt ...


      // change trace start time of stations
      for (auto sIt = rEvent.StationsBegin(); sIt != rEvent.StationsEnd(); ++sIt) {
        Station& station = *sIt;

        // apply beacon correction only if time offset could be determined successfully
        if (timeOffsetUncertainty[station.GetId()] < 99.) { // is set to 99 correction fails
          // station.SetTraceStartTime(station.GetTraceStartTime() + timeOffsets[station.GetId()]);
          station.GetRecData().SetParameter(eTraceStartTime, station.GetRecData().GetParameter(eTraceStartTime) + timeOffsets[station.GetId()], false);
          station.GetRecData().SetParameterError(eTraceStartTime,fBeaconTimeCalibrationUncertainty, false);
          // store value that has been added to the trace start time
          station.GetRecData().SetParameter(eTraceStartTimeCorrectionByBeacon, timeOffsets[station.GetId()]);
          // also change relative position of signal search window accordingly
          station.GetRecData().SetParameter(eSignalSearchWindowStart, station.GetRecData().GetParameter(eSignalSearchWindowStart) - timeOffsets[station.GetId()], false);
          // also change relative position of signal search window accordingly
          station.GetRecData().SetParameter(eSignalSearchWindowStop, station.GetRecData().GetParameter(eSignalSearchWindowStop) - timeOffsets[station.GetId()], false);

          message.str("");
          message << "Shifted trace starttime of station "
                  << station.GetId()
                  << " by "
                  << timeOffsets[station.GetId()]
                  << " ns; \tuncert. ~ "
                  << timeOffsetUncertainty[station.GetId()]
                  << " ns";
          if (fInfoLevel >= eInfoIntermediate)
            INFO(message.str());

          // add statistics
          --fBeaconCorrectionFailures[station.GetId()]; // at the beginning a failure is assumed, and substracted, if this point is reached
          ++fBeaconCorrectionSuccesses[station.GetId()];
        } else {
          if (fRejectIfCorrectionFailed)
            station.SetRejectedReason(eBeaconCorrectionFailed);
          message.str("");
          message << "Shifted trace starttime of station "
                  << station.GetId()
                  << " by 0 ns - Beacon correction failed!";
                  if (fRejectIfCorrectionFailed)
                    message << " Station will be rejected.\n";
                  else
                    message << "\n.";
          if (fInfoLevel >= eInfoIntermediate)
            INFO(message.str());
        }
      }

      // equalize time stamps of stations (if different, and if wanted)
      // by multiplying a phase gradient to each spectrum
      if (fEqualizeStationTimeStamps)
        matchStationTimeStamps(rEvent);

      // write debug output to file, if filename is given in XML file
      if (fDebugStationIds.size() > 0)
        if (!writeDebugFile(refStationID, timeOffsets, timeOffsetUncertainty, rEvent.GetHeader().GetTime()+ refStationTimeStamp)) {
          message.str("");
          message << "Error while writing debug file: "
                  << fDebugOutputFile;
          ERROR(message.str());
          return eFailure;
        }

    }

    // remember that module already has been called (important for debug output)
    fFirstModuleCall = false;

    return eSuccess;
  }


  VModule::ResultFlag
  RdChannelBeaconTimingCalibrator::Finish()
  {
    INFOFinal("Finish()");

    ostringstream info;
    info << "Statistics of Beacon correction:\n"
            "Station ID       #successes    #failures   SNR(freq1)   SNR(freq2)   SNR(freq3)   SNR(freq4)\n"
            "--------------------------------------------------------------------------------------------\n";

    for (const auto& ie : fBeaconCorrectionFailures) {
      const auto id = ie.first;
      info << "  " << id << "          \t"
           << fBeaconCorrectionSuccesses[id] << "  \t      "
           << fBeaconCorrectionFailures[id] << "  \t  "
           << fBeaconSNRs[id][0] / fBeaconSNRcounter[id] << "  \t  "
           << fBeaconSNRs[id][1] / fBeaconSNRcounter[id] << "  \t  "
           << fBeaconSNRs[id][2] / fBeaconSNRcounter[id] << "  \t  "
           << fBeaconSNRs[id][3] / fBeaconSNRcounter[id] << '\n';
    }

    info << "--------------------------------------------------------------------------------------------";
    INFO(info);
    return eSuccess;
  }


  std::vector<std::map<int, double> >
  RdChannelBeaconTimingCalibrator::getReferencePhases(std::vector<double> beaconFrequencies)
  {
    static std::vector<std::map<int, double> > referencePhases;

    // check, if reference phases have already been read from file
    if (referencePhases.size() > 0)
      return referencePhases;

    // check if name of text file with reference phases is provided,
    // otherwise read reference phase from database

    // Create file name containing frequency
    ifstream phasefile;
    phasefile.open(fInputPhaseFile.c_str());

    stringstream message;
    if ((phasefile.is_open()) && (phasefile.good())) {
      message.str("");
      message << "Reading reference phases from file: "
              << fInputPhaseFile;
      if (fInfoLevel >= eInfoDebug)
        INFO(message.str());
    } else {
      message.str("");
      message << "Cannot open file for reference phases: "
              << fInputPhaseFile;
      ERROR(message.str());
      return referencePhases;
    }

    // read file header to check consistency
    stringstream headerDummy;
    string readString;
    phasefile >> readString;

    if (readString != "Station") {
      ERROR("Wrong format of file header");
      phasefile.close();
      return referencePhases;
    }

    for (unsigned int beaconFreqN = 0; beaconFreqN < beaconFrequencies.size(); ++beaconFreqN) {
      if (!phasefile.good() ) {
        ERROR("Wrong format of file header");
        return referencePhases;
      }
      // test header content (do beacon frequencies match?)
      phasefile >> readString;
      headerDummy.str("");
      headerDummy << "RefPhase[rad]_for_" << beaconFrequencies[beaconFreqN]/megahertz << "MHz";
      if (readString != headerDummy.str()) {
        if (fBeaconFrequencies.size() < 1)
          ERROR("Beacon frequencies in file header do not match frequencies read from database! Maybe only the order of the frequencies is wrong.");
        else
          ERROR("Beacon frequencies in file header do not match XML file!");
        phasefile.close();
        return referencePhases;
      }
    }

    // loop through all lines and read the station and the reference phases
    referencePhases = vector<map<int, double> > (beaconFrequencies.size()); // delete previous content of vector

    while (phasefile.good()) {
      unsigned int stationID=0;
      phasefile >> stationID;
      // check if there is still a line in the file
      if (!phasefile.good())
        break;
      for (unsigned int beaconFreqN = 0; beaconFreqN < beaconFrequencies.size(); ++beaconFreqN) {
        double refPhase = 0;
        phasefile >> refPhase;
        referencePhases[beaconFreqN][stationID] = refPhase;
      }
    }

    phasefile.close();

    return referencePhases;
  }


  void
  RdChannelBeaconTimingCalibrator::matchStationTimeStamps(revt::REvent& rEvent) const
  {
    // store time stamp of the first station as reference
    double firstTimeStamp  = rEvent.StationsBegin()->GetRecData().GetParameter(eTraceStartTime);
    // loop through all stations:
    // shift traces if the time stamp does not match the one of the first station
    for (auto sIt = rEvent.StationsBegin(); sIt != rEvent.StationsEnd(); ++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);
        }

        // set trace start time, and change relative position of signal search window accordingly
        // since the signal and the trace start time are shifted consistingly, there is no correction value stored for the trace start time
        station.GetRecData().SetParameter(eTraceStartTime,firstTimeStamp, false);
        station.GetRecData().SetParameter(eSignalSearchWindowStart, station.GetRecData().GetParameter(eSignalSearchWindowStart) + (stationTimeStamp-firstTimeStamp), false);
        station.GetRecData().SetParameter(eSignalSearchWindowStop, station.GetRecData().GetParameter(eSignalSearchWindowStop) + (stationTimeStamp-firstTimeStamp), false);
      }
    }
  }


  bool
  RdChannelBeaconTimingCalibrator::writeDebugFile(const int& refStationID,
                                                  const std::map<int,double>& timeOffsets,
                                                  const std::map<int,double>& timeUncertainties,
                                                  const utl::TimeStamp& eventGPStimeStamp) const
  {
    // When the beacon module is called for the first time, the file should not exist.
    // Thus, test if file exists, and if not, create file header
    if (fFirstModuleCall) {
      if (!fOverwriteDebugFile) {
        ifstream testReading;
        testReading.open(fDebugOutputFile.c_str());
        if (testReading.is_open()) {
          stringstream message;
          message << "Debug file '"
                  << fDebugOutputFile
                  << "' already exists! No debug output written!\n";
          ERROR(message.str());
          return false; // found error
        }
        testReading.close();
      }

      // create file header
      ofstream debugFileHeaderWriting(fDebugOutputFile.c_str());
      debugFileHeaderWriting << "GPStime[s]          \tRefStationID";
      for (set<int>::iterator it=fDebugStationIds.begin(); it != fDebugStationIds.end(); ++it)
        debugFileHeaderWriting << " \t" << *it << "_offset[ns] \t" << *it << "_uncert[ns]";
      debugFileHeaderWriting << endl;
      debugFileHeaderWriting.close();
    }

    // append event to file
    ofstream outputfile;
    outputfile.open(fDebugOutputFile.c_str(), ofstream::app);

    // check if file could be opened
    if (!(outputfile.is_open())) {
      stringstream message;
      message << "Failed to open beacon debug file: "
              << fDebugOutputFile;
      ERROR(message.str());
      return false;
    }

    // time stamp and reference station ID
    outputfile << eventGPStimeStamp.GetGPSSecond() << "."
               << setfill ('0') << setw(9) << int(round(eventGPStimeStamp.GetGPSNanoSecond()))
               << " \t" << refStationID;

    // loop through station IDs and write offsets and uncertainties to file
    for (auto it = fDebugStationIds.begin(); it != fDebugStationIds.end(); ++it) {
      // check if station exists for this event
      if ( timeUncertainties.find(*it) != timeUncertainties.end() ) {
        // check if beacon correction worked
        if ( timeUncertainties.find(*it)->second != 99. )
          outputfile << " \t" << setfill (' ') << setw(7) << timeOffsets.find(*it)->second
                     << " \t" << setfill (' ') << setw(7) << timeUncertainties.find(*it)->second;
        else
          outputfile << " \t" << "  xxx  " << " \t" << "  xxx  ";
      } else {
        outputfile << " \t" << "  ---  " << " \t" << "  ---  ";
      }
    }
    outputfile << endl;
    outputfile.close();

    return true; // everything ok
  }

}