FindTimeOffset.cc

Go to the documentation of this file.

#include <evt/Event.h>
#include <evt/ShowerFRecData.h>

#include <fdet/FDetector.h>
#include <fdet/Eye.h>
#include <fevt/FEvent.h>
#include <fevt/Eye.h>
#include <fevt/EyeHeader.h>
#include <fevt/EyeRecData.h>

#include <utl/ErrorLogger.h>
#include <utl/PhysicalConstants.h>
#include <utl/TimeStamp.h>
#include <utl/Transformation.h>
#include <utl/GeometryUtilities.h>
#include <utl/FileName.h>

#include <fwk/LocalCoordinateSystem.h>
#include <fwk/CoordinateSystemRegistry.h>

#include <sdet/SDetector.h>
#include <sdet/Station.h>
#include <sevt/SEvent.h>
#include <sevt/Station.h>

#include <sstream>
#include <fstream>

#include "FindTimeOffset.h"


namespace FindTimeOffset {

  fwk::VModule::ResultFlag
  FindTimeOffset::Run(evt::Event& event)
  {
    if (!event.HasFEvent())
      return fwk::VModule::eSuccess;

    if (!event.HasSEvent()) {
      INFO("No SEvent");
      return fwk::VModule::eContinueLoop;
    }
    auto& sEvent = event.GetSEvent();

    bool isCeleste = false;
    bool isRamiro = false;

    for (const auto& station : sEvent.StationsRange()) {
      if (!station.HasRecData() || !station.IsCandidate())
        continue;
      if (station.GetId() == 203)
        isCeleste = true;
      else if (station.GetId() == 805)
        isRamiro = true;
    }

    if (isCeleste)
      INFO("This is a CLF event.");
    if (isRamiro)
      INFO("This is an XLF event.");
    if (isCeleste && isRamiro) {
      INFO("Cannot have both Celeste and Ramiro; skipping.");
      return fwk::VModule::eContinueLoop;
    }
    if (!isCeleste && !isRamiro) {
      INFO("No Celeste or Ramiro found: NOT a hybrid laser event; skipping.");
      return fwk::VModule::eContinueLoop;
    }

    ExamineFEvent(event, isCeleste);

    return fwk::VModule::eSuccess;
  }


  void
  FindTimeOffset::ExamineFEvent(const evt::Event& event, const bool isCeleste)
  {
    const auto& sEvent = event.GetSEvent();
    const auto& fEvent = event.GetFEvent();
    for (auto eye = fEvent.EyesBegin(fevt::ComponentSelector::eHasData);
         eye != fEvent.EyesEnd(fevt::ComponentSelector::eHasData); ++eye)
      FindOffset(sEvent, *eye, isCeleste);
  }


  template<class V>
  std::string
  ToString(const V& v, const utl::CoordinateSystemPtr& cs, const double unit)
  {
    std::ostringstream os;
    os << '('
       << v.GetX(cs)/unit << ", "
       << v.GetY(cs)/unit << ", "
       << v.GetZ(cs)/unit << ')';
    return os.str();
  }


  void
  FindTimeOffset::FindOffset(const sevt::SEvent& sEvent, const fevt::Eye& eye,
                             const bool isCeleste)
  {
    const utl::CoordinateSystemPtr& laserCs = (isCeleste) ?
      fwk::CoordinateSystemRegistry::Get("CLF") :
      fwk::CoordinateSystemRegistry::Get("XLF");
    const utl::Point laserPos(0,0,0, laserCs);

    const auto& det = det::Detector::GetInstance();
    const auto& dEye = det.GetFDetector().GetEye(eye);
    const auto& eyePos = dEye.GetPosition();
    const auto& eyeCs = dEye.GetLocalCoordinateSystem();
    const utl::Vector verticalAtEye(0,0,1, eyeCs);
    const auto eye2Laser = laserPos - eyePos;

    const auto& eyeRec = eye.GetRecData();

    const double chi0 = eyeRec.GetChiZero();
    const double rp = eyeRec.GetRp();
    const double t0 = eyeRec.GetTZero();
    const auto& sdp = eyeRec.GetSDP();

    const auto& fdEventTime = eye.GetHeader().GetTimeStamp();
    // DV: we should look this up from fdet
    const utl::TimeInterval fdTraceLength = 1000 * 100*utl::ns;
    // DV: absolute reference for all relative times (ie when given as double,
    // eg t0 or centroids etc) is the EyeHeader::TimeStamp minus the trace length
    const auto fdEyeTriggerTime = fdEventTime - fdTraceLength;

    // DV: moved the *= -1 in original code to swap the cross product
    const utl::Vector sdpHorizontalAtEye = Normalized(Cross(verticalAtEye, sdp));

    // obtain laser pointing by rotating horizontal Vector around the sdp normal by -chi0
    const auto laserPointing =
      utl::Transformation::Rotation(-chi0, sdp, eyeCs) * sdpHorizontalAtEye;

    std::ostringstream info;
    info << "GPS = " << fdEventTime.GetGPSSecond() << " s "
         << int(fdEventTime.GetGPSNanoSecond()) << " ns, "
            "eye = " << eye.GetId() << ", "
            "chi0 = " << chi0/utl::deg << " deg, "
            "Rp = " << rp/utl::meter << " m, "
            "t0 = " << t0/utl::ns << " ns";
    INFO(info);

    info.str("");
    info << "laser pointing (eye CS): "
            "theta = " << laserPointing.GetTheta(eyeCs)/utl::deg << " deg, "
            "phi = " << laserPointing.GetPhi(eyeCs)/utl::deg << " deg";
    INFO(info);

    info.str("");
    info << "laser pointing (laser CS): "
            "theta = " << laserPointing.GetTheta(laserCs)/utl::deg << ", "
            "phi = " << laserPointing.GetPhi(laserCs)/utl::deg << " deg";
    INFO(info);

    const double beta = utl::Angle(-eye2Laser, laserPointing);

    // DV: delays have to be subtracted from time
    // geometrical delay is negative since it's early
    const utl::TimeInterval geometricDelay = -eye2Laser.GetMag() * std::cos(beta) / utl::kSpeedOfLight;
    const utl::TimeInterval fiberDelay = 180*utl::ns;

    info.str("");
    info << "beta = " << beta/utl::deg << " deg, "
            "dt = " << geometricDelay/utl::ns << " ns";
    INFO(info);

    const auto& shFRec = eyeRec.GetFRecShower();
    const auto& fdCore = shFRec.GetCorePosition();
    const auto& fdCoreTime = shFRec.GetCoreTime();
    const auto& fdAxis = shFRec.GetAxis();

    for (const auto& station : sEvent.StationsRange()) {

      if (!station.HasRecData() ||
          !station.IsCandidate() ||
          !(station.GetId() == 203 || station.GetId() == 805))
        continue;

      // the offset is obtained as the difference of sd and fd times in the t0 point
      const auto& stationStart = station.GetRecData().GetSignalStartTime();
      const auto laserStart = stationStart - fiberDelay;
      // time at the t0 point from sd
      const auto sdT0 = laserStart - geometricDelay;

      info.str("");
      info << "sdT0 = " << sdT0;
      INFO(info);

      // time at the t0 point from fd
      const auto fdT0 = fdEyeTriggerTime + t0;

      info.str("");
      info << "fdT0 = " << fdT0;
      INFO(info);

      const auto offset = sdT0 - fdT0;

      info.str("");
      info << "SD/FD time offset = " << offset/utl::ns << " ns";
      INFO(info);

      // the alternative offset2 is the difference of the sd time and the time
      // of the fd-reconstructed shower front arriving at the xlf/clf

      // correct the core time by a (small) axis-parallel shift delay for when the front passes through the station
      const double coreLaserParallelDistance = fdAxis * (fdCore - laserPos);
      const double coreLaserParallelDt = coreLaserParallelDistance / utl::kSpeedOfLight;
      const auto fdTimeAtLaser = fdCoreTime + coreLaserParallelDt;

      const auto offset2 = laserStart - fdTimeAtLaser;

      const auto& stationPos = det.GetSDetector().GetStation(station).GetPosition();
      const auto projectedFdCore =
        utl::Intersection(utl::Line(fdCore, fdAxis), utl::Plane(laserPos, utl::Vector(0,0,1, laserCs)));
      info.str("");
      info << "(in CS of the laser)\n"
              "  fdCore = " << ToString(fdCore, laserCs, utl::m) << " m\n"
              "  fdCore_proj = " << ToString(projectedFdCore, laserCs, utl::m) << " m\n"
              "  stationPos = " << ToString(stationPos, laserCs, utl::m) << " m\n"
              "  laserPos = " << ToString(laserPos, laserCs, utl::m) << " m\n"
              "  distance(core, laser) = " << (fdCore - laserPos).GetMag()/utl::meter << " m\n"
              "  distance(core, station) = " << (fdCore - stationPos).GetMag()/utl::meter << " m\n"
              "  laserStart - fdCoreTime = " << (laserStart - fdCoreTime)/utl::ns << " ns\n"
              "  projection delay = " << coreLaserParallelDt/utl::ns << " ns";
      INFO(info);

      info.str("");
      info << "SD/FD time offset from core = " << offset2/utl::ns << " ns";
      INFO(info);

      const double laserPointingElevation = utl::kHalfPi - laserPointing.GetTheta(laserCs);

      const auto outputFile = "offset.dat";
      if (!utl::File::Exists(outputFile)) {
        std::ofstream fout(outputFile);
        fout << "# fd_id offset gps_sec gps_ns laser_elevation laser_azimuth signal celeste offset2 fd_core_x y z\n";
      }
      std::ofstream fout(outputFile, std::ios::app);
      fout << eye.GetId() << ' '
           << offset/utl::ns << ' '
           << fdEventTime.GetGPSSecond() << ' '
           << int(fdEventTime.GetGPSNanoSecond()) << ' '
           << laserPointingElevation/utl::deg << ' '
           << laserPointing.GetPhi(laserCs)/utl::deg << ' '
           << station.GetRecData().GetTotalSignal() << ' '
           << isCeleste << ' '
           << offset2/utl::ns << ' '
           << projectedFdCore.GetX(laserCs)/utl::m << ' '
           << projectedFdCore.GetY(laserCs)/utl::m << ' '
           << projectedFdCore.GetZ(laserCs)/utl::m
           << '\n';

    }
  }

}