SdStationPositionCorrection.cc

Go to the documentation of this file.

#include "SdStationPositionCorrection.h"
#include <evt/Event.h>
#include <sevt/SEvent.h>
#include <sdet/SDetector.h>
#include <fwk/LocalCoordinateSystem.h>
#include <fwk/CentralConfig.h>
#include <utl/ReadStream.h>
#include <utl/UTMPoint.h>
#include <utl/LeapSeconds.h>
#include <utl/Math.h>
#include <utl/zfstream.h>
#include <utl/String.h>

using namespace evt;
using namespace sevt;
using namespace std;
using namespace utl;
using namespace fwk;


namespace SdStationPositionCorrectionOG {

  class UnrestrictedSDetStation : public sdet::Station {
  public:
    void
    SetPosition(const Point& pos,
                const CoordinateSystemPtr& referenceSystem = CoordinateSystemPtr())
    {
      fPosition = pos;
      fReferenceSystem = referenceSystem;
    }
  };


  struct UTM {

    double fNorthing;
    double fEasting;
    double fHeight;

    friend istream& operator>>(istream& is, UTM& u)
    { return is >> u.fNorthing >> ws >> u.fEasting >> ws >> u.fHeight; }

  };


  struct KnownPositions {

    unsigned int fStationId;
    UTM fPMS;
    UTM fGPSAverage;
    UTM fGPSAverageWithDOPWeights;

    friend istream& operator>>(istream& is, KnownPositions& p)
    { return is >> p.fStationId >> ws >> p.fPMS >> ws >> p.fGPSAverage >> ws >> p.fGPSAverageWithDOPWeights; }

  };


  struct FixedPosition {

    unsigned short int fStationId;
    time_t fUnixSecond;
    int fNRepeats;
    UTM fPosition;

    friend istream& operator>>(istream& is, FixedPosition& p)
    { return is >> p.fStationId >> ws >> p.fUnixSecond >> ws >> p.fNRepeats >> ws >> p.fPosition; }

  };


  VModule::ResultFlag
  SdStationPositionCorrection::Init()
  {
    // trigger detector update to read standard station positions
    auto& det = det::Detector::GetInstance();
    det.Update(GetCurrentSystemTime());
    const auto& sDet = det.GetSDetector();

    auto topB = CentralConfig::GetInstance()->GetTopBranch("SdStationPositionCorrection");

    fEffectiveVerticalDelay = topB.GetChild("EffectiveVerticalDelay").Get<double>();

    const auto var = topB.GetChild("EffectiveVarianceCoefficients").Get<vector<double>>();
    fEffectiveVarianceCoefficients = Triple(var.at(0), var.at(1), var.at(2));

    const auto sdStationGPSPositionsFile = topB.GetChild("SdStationGPSPositionsFile").Get<string>();

    const auto sdStationFixPositionsFile = topB.GetChild("SdStationFixPositionsFile").Get<string>();

    const int zone = 19;
    const char band = 'H';
    const auto& ellipsoid = ReferenceEllipsoid::GetWGS84();
    const auto& siteCS = det.GetSiteCoordinateSystem();

    typedef std::map<unsigned int, utl::Point> PositionMap;
    PositionMap gpsPositions;

    {
      zifstream posFile(sdStationGPSPositionsFile.c_str());

      if (!posFile.is_open()) {
        ostringstream err;
        err << "Cannot open GPS positions file '" << sdStationGPSPositionsFile << "'!";
        ERROR(err);
        return eFailure;
      }

      vector<KnownPositions> positions;
      ReadStream(posFile).GetLines(positions, /*skip comments*/true, /*trim*/true);

      ostringstream info;
      info << "Got " << positions.size() << " GPS positions from file \"" << sdStationGPSPositionsFile << "\"";
      INFO(info);

      for (const auto& kpos : positions) {
        const unsigned int sId = kpos.fStationId;
        try {
          const auto& dStation = sDet.GetAllStation(sId);
          const auto& detPosition = dStation.GetPosition();
          const auto& gps = kpos.fGPSAverage;
          // skip entries with only PMS positions
          if (gps.fNorthing && gps.fEasting && gps.fHeight) {
            const auto position = UTMPoint(gps.fNorthing, gps.fEasting, gps.fHeight, zone, band, ellipsoid).GetPoint(siteCS);
            gpsPositions.insert(make_pair(sId, position));
            const auto& localCS = LocalCoordinateSystem::Create(position);
            static_cast<UnrestrictedSDetStation&>(const_cast<sdet::Station&>(dStation)).SetPosition(position, localCS);
          }
        } catch (NonExistentComponentException&) {
          ostringstream warn;
          warn << "Station id=" << sId << " has a GPS position but is not found in the SDetector!";
          WARNING(warn);
        }
      }
    }

    fStationSpatialOffset.clear();

    {
      zifstream fixFile(sdStationFixPositionsFile.c_str());

      if (!fixFile.is_open()) {
        ostringstream err;
        err << "Cannot open fix-positions file '" << sdStationFixPositionsFile << "'!";
        ERROR(err);
        return eFailure;
      }

      vector<FixedPosition> fixPositions;
      ReadStream(fixFile).GetLines(fixPositions, /*strip comments*/true, /*trim*/true);

      ostringstream info;
      info << "Got " << fixPositions.size() << " fixed positions from file \"" << sdStationFixPositionsFile << "\"";
      INFO(info);

      set<int> missingStations;
      for (const auto& fpos : fixPositions) {
        const unsigned int sId = fpos.fStationId;
        auto& st = fStationSpatialOffset[sId];
        unsigned long gpsSec = 0;
        LeapSeconds::GetInstance().ConvertUnixToGPS(fpos.fUnixSecond, gpsSec);
        auto& dx = st[gpsSec];
        dx = Triple(0, 0, 0);
        // gps measurement campaings will have negative NRepeats
        if (fpos.fNRepeats > 0) { // for gps measurement campaigns assume correlated changes in time
          const auto gpsIt = gpsPositions.find(sId);
          if (gpsIt == gpsPositions.end())
            missingStations.insert(sId);
          else {
            const auto& fix = fpos.fPosition;
            const auto pos = UTMPoint(fix.fNorthing, fix.fEasting, fix.fHeight, zone, band, ellipsoid).GetPoint(siteCS);
            const auto& gpsPos = gpsIt->second;
            const auto diff = gpsPos - pos;
            const auto& localCS = LocalCoordinateSystem::Create(gpsPos);
            dx = diff.GetCoordinates(localCS);
          }
        }
      }

      if (!missingStations.empty()) {
        ostringstream warn;
        warn << "Stations without GPS positions will have delay set to zero: " << String::OfSortedIds(missingStations);
        WARNING(warn);
      }
    }

    return eSuccess;
  }


  VModule::ResultFlag
  SdStationPositionCorrection::Run(Event& event)
  {
    if (!event.HasSEvent())
      return eSuccess;

    auto& sEvent = event.GetSEvent();

    vector<pair<unsigned int, pair<double, double>>> stationDtSigma;
    for (auto& station : sEvent.StationsRange()) {
      if (!station.HasRecData() || station.IsUUB())
        continue;
      auto& recData = station.GetRecData();
      const auto& oldTime = recData.GetSignalStartTime();
      if (!oldTime)
        continue;
      const unsigned int sId = station.GetId();
      const auto tdxIt = fStationSpatialOffset.find(sId);
      if (tdxIt == fStationSpatialOffset.end())
        continue;
      const auto& tdx = tdxIt->second;
      if (tdx.empty())
        continue;
      auto dxIt = tdx.lower_bound(oldTime.GetGPSSecond());
      if (dxIt == tdx.begin())
        continue;
      --dxIt;
      const auto& dx = dxIt->second;
      const TimeInterval dt = dx.get<2>() * fEffectiveVerticalDelay;  // Z
      recData.SetStartTimeGPSCorrection(dt);
      recData.SetSignalStartTime(oldTime - dt);
      const double sigma2 =
        Sqr(fEffectiveVarianceCoefficients.get<0>() * dx.get<0>()) +  // X
        Sqr(fEffectiveVarianceCoefficients.get<1>() * dx.get<1>()) +  // Y
        Sqr(fEffectiveVarianceCoefficients.get<2>() * dx.get<2>());   // Z
      recData.SetGPSTimeVariance(sigma2);
      stationDtSigma.push_back(make_pair(sId, make_pair(dt, sqrt(sigma2))));
    }

    if (stationDtSigma.empty())
      INFO("no stations were corrected.");
    else {
      ostringstream info;
      info << "station time offset and sigma [ns]:";
      for (const auto& ids : stationDtSigma)
        info << " (" << ids.first << ": "
             << ids.second.first/nanosecond << " +- " << ids.second.second/nanosecond << ')';
      INFO(info);
    }

    return eSuccess;
  }

}