FOVCalculator.cc

Go to the documentation of this file.

#include <utl/config.h>
#include <adst/RecEvent.h>

#include "FOVCalculator.h"

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

#include <utl/Reader.h>
#include <utl/TimeStamp.h>
#include <utl/ErrorLogger.h>
#include <utl/MathConstants.h>
#include <utl/PhysicalConstants.h>
#include <utl/Transformation.h>
#include <utl/CoordinateSystemPtr.h>
#include <utl/Point.h>
#include <utl/Vector.h>
#include <utl/AxialVector.h>
#include <utl/UTMPoint.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/TabulatedFunction.h>

#include <det/Detector.h>
#include <fdet/FDetector.h>
#include <fdet/Pixel.h>
#include <fdet/Camera.h>
#include <fdet/Eye.h>
#include <fdet/Telescope.h>

#include <fevt/FEvent.h>
#include <fevt/Eye.h>
#include <fevt/Telescope.h>

#include <atm/ProfileResult.h>
#include <atm/InclinedAtmosphericProfile.h>
#include <atm/AttenuationResult.h>

#include <iostream>
#include <string>
#include <map>
#include <vector>

using namespace fwk;
using namespace utl;
using namespace fdet;
using namespace atm;
using namespace std;
using namespace otoa;


void
FOVCalculator::FillFOVVariables(const fevt::FEvent& fdEvent, RecEvent& theRecEvent)
{
  CoordinateSystemPtr referenceCS = det::Detector::GetInstance().GetReferenceCoordinateSystem();
  const Atmosphere& atm = det::Detector::GetInstance().GetAtmosphere();
  const fdet::FDetector& fdDet  = det::Detector::GetInstance().GetFDetector();

  std::vector<FDEvent>& theFDEvents = theRecEvent.GetFDEvents();
  for (unsigned int i = 0; i < theFDEvents.size(); ++i) {

    const unsigned int eyeID = theFDEvents[i].GetEyeId();
    const fdet::Eye& detEye = fdDet.GetEye(eyeID);
    const fevt::Eye& evtEye = fdEvent.GetEye(eyeID, fevt::ComponentSelector::eInDAQ);

    set<unsigned int> telescopesWithData;
    set<unsigned int> telescopesInDAQ;
    for (auto telIter = evtEye.TelescopesBegin(fevt::ComponentSelector::eInDAQ);
         telIter != evtEye.TelescopesEnd(fevt::ComponentSelector::eInDAQ); ++telIter)
      telescopesInDAQ.insert(telIter->GetId());

    for (auto telIter = evtEye.TelescopesBegin(fevt::ComponentSelector::eHasData);
         telIter != evtEye.TelescopesEnd(fevt::ComponentSelector::eHasData); ++telIter)
      telescopesWithData.insert(telIter->GetId());

    if (theFDEvents[i].GetRecLevel() >= eHasGHParameters) {

      FdRecShower& theShower = theFDEvents[i].GetFdRecShower();
      const double zeta = theFDEvents[i].GetFdRecApertureLight().GetZeta();
      const TVector3& adstCore = theShower.GetCoreSiteCS();
      const TVector3& adstAxis = theShower.GetAxisSiteCS();

      const Vector axis(adstAxis.X(), adstAxis.Y(), adstAxis.Z(), referenceCS);
      const Point core(adstCore.X(), adstCore.Y(), adstCore.Z(), referenceCS);

      try {

        atm.InitSlantProfileModel(core, axis, 50*g/cm2);

        double xLow = -1;
        double xUp = -1;

        CalculateFOVBoundaries(core, axis, detEye,
                               telescopesInDAQ,
                               telescopesWithData,
                               zeta, xLow, xUp);

        theShower.SetXFOVMin(xLow);
        theShower.SetXFOVMax(xUp);

      } catch (InclinedAtmosphericProfile::InclinedAtmosphereModelException& e) {
        ERROR(e.GetMessage());
      }
    }
  }

  CalculateXmaxErrors(theRecEvent);
}


void
FOVCalculator::CalculateXmaxErrors(RecEvent& theRecEvent)
{
  fXmaxScanner.EstimateXmaxErrors(theRecEvent);
  std::vector<FDEvent>& theFDEvents = theRecEvent.GetFDEvents();
  const std::vector<LongitudinalScan> longScan =
    fXmaxScanner.GetLongitudinalScan();
  const std::vector<double> errAtXmaxVec =
    fXmaxScanner.GetErrorAtXmax();
  if (theFDEvents.size() != longScan.size() ||
      theFDEvents.size() != errAtXmaxVec.size()) {
    ERROR("size mismatch FDEvent vs. XmaxErrors!!");
    cout << theFDEvents.size() << ' '
         << longScan.size() << ' '
         << errAtXmaxVec.size() << endl;
    throw AugerException();
  }

  for (unsigned int i = 0; i < theFDEvents.size(); ++i) {
    FdRecShower& recShower = theFDEvents[i].GetFdRecShower();
    vector<double>& zVec = recShower.GetZVector();
    zVec.clear();
    vector<double>& mvaVec = recShower.GetMinViewingAngles();
    mvaVec.clear();
    const vector<double>& xmaxErrors = longScan[i].GetXiXmaxVector();
    const vector<double>& mvaValues = longScan[i].GetMinViewingAngleVector();
    for (unsigned int j = 0; j < xmaxErrors.size(); ++j) {
      zVec.push_back(xmaxErrors[j]/g*cm2);
      mvaVec.push_back(mvaValues[j]);
    }
    recShower.SetZXmax(errAtXmaxVec[i]/g*cm2);
  }
}


bool
FOVCalculator::CalculateFOVBoundaries(const Point& core, const Vector& axis,
                                      const fdet::Eye& detEye,
                                      const set<unsigned int>& telsInDAQ,
                                      const set<unsigned int>& telsWithData,
                                      const double zeta,
                                      double& xLow, double& xUp)
{
  const Atmosphere& atm = det::Detector::GetInstance().GetAtmosphere();
  ReferenceEllipsoid wgs84(ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84));

  const Point& eyePos = detEye.GetPosition();
  CoordinateSystemPtr eyeCS = detEye.GetEyeCoordinateSystem();

  double xUpGeom = -1;
  double xLowGeom = -1;
  double xUpDAQ = -1;
  double xLowDAQ = -1;

  try {

    const atm::ProfileResult& depthProfile = atm.EvaluateSlantDepthVsDistance();
    const double maxDistance = depthProfile.MaxX();
    const double minDistance = depthProfile.MinX();
    const CoordinateSystemPtr& CS = detEye.GetLocalCoordinateSystem();
    const utl::Vector vertical(0, 0, 1, CS, Vector::kCartesian);

    const double cosZeta = cos(zeta);

    const utl::Vector eyeToShower = core - eyePos;

    utl::Vector sdp = Cross(axis, eyeToShower);
    const double Rp = sdp.GetMag();
    sdp.Normalize();
    sdp.TransformTo(eyeCS);

    const AxialVector horizontal = Normalized(Cross(sdp, vertical));

    const double chi0 = Angle(horizontal, axis);
    const Vector eyeCoreVec = Rp / sin(kPi - chi0) * horizontal;
    const double eyeCoreDistance = eyeCoreVec.GetMag();

    for (auto telIter = detEye.TelescopesBegin(); telIter != detEye.TelescopesEnd(); ++telIter) {

      const unsigned int telID = telIter->GetId();

      const bool telHasData = (telsWithData.find(telID) != telsWithData.end());
      const bool telInDAQ = (telsInDAQ.find(telID)!= telsInDAQ.end());

      bool goodTelescope = false;

      try {
        const bool telInAquisition = telInDAQ;
        const double telUpTime = telIter->GetUpTimeFraction();
        if (fVerbosity > -1)
          cout << " telId=" << telID << " upTime = " << telUpTime
               << " isInAqu=" << telInAquisition
               << " " << (telHasData ? "[data]" : "      ")
               << '\n';
        goodTelescope = (telUpTime > 0 && telInDAQ);
      } catch (NonExistentComponentException& c) {
        // use telescopes in event for non-existing uptime
        cerr << " no FdUptime available!!\n"
             << c.GetExceptionName() << "\n " << c.GetMessage() << '\n';
        if (telHasData)
          goodTelescope = true;
      }

      if (!goodTelescope && telHasData) {
        cerr << " FOVCalculator::CalculateFOVBoundaries() -"
             << " telescope out of DAQ with data!! eye="
             << detEye.GetId() << " tel=" << telIter->GetId() << '\n';
        goodTelescope = true;
      }

      if (!goodTelescope)
        continue;

      vector<unsigned int> candidatePixelIds;

      for (unsigned int iP = telIter->GetFirstPixelId(); iP <= telIter->GetLastPixelId(); ++iP) {

        const Vector& pixelDir = telIter->GetPixel(iP).GetDirection();

        const double angle = fabs(0.5*kPi - Angle(sdp, pixelDir));
        if (angle <= 0.75*degree) {
          candidatePixelIds.push_back(iP);
        }
      }

      const unsigned int minSLTPixels = 4;

      if (candidatePixelIds.size() >= minSLTPixels) {

        for (unsigned int iPix = 0; iPix < candidatePixelIds.size(); ++iPix) {

          const unsigned int iP = candidatePixelIds[iPix];
          const Vector& pixelDir = telIter->GetPixel(iP).GetDirection();

          const Vector chi_i_vector = pixelDir - (sdp * pixelDir) * sdp;

          const double nominalAngle = Angle(chi_i_vector, horizontal);
          if (nominalAngle < 0 || nominalAngle>chi0)
            continue;

          for (double dChi = -0.85*degree; dChi <= 0.85*degree; dChi += 0.1*degree) {

            const double chi_i = nominalAngle + dChi;

            double l = eyeCoreDistance * sin(chi_i) / sin(chi0 - chi_i);
            if (l < 0)
              l *= -1;
            double distance = -l; // by convention
            if (distance < minDistance)
              distance = minDistance;
            else if (distance > maxDistance)
              distance = maxDistance;

            const Vector closestApproach = eyeCoreVec + l * axis;
            closestApproach.TransformTo(eyeCS);

            // check if close to border (FdApertureLightFinder::IsNearBorder())
            bool isNearBorder = false;

            for (auto oobPixIter = telIter->OutOfBorderPixelsBegin(); oobPixIter!= telIter->OutOfBorderPixelsEnd(); ++oobPixIter) {
              if (CosAngle(closestApproach, *oobPixIter) > cosZeta) {
                isNearBorder = true;
                break;
              }
            }

            if (!isNearBorder) {

              const Point pclosestApproach(closestApproach.GetX(eyeCS),
                                           closestApproach.GetY(eyeCS),
                                           closestApproach.GetZ(eyeCS),
                                           eyeCS);
              const double thePhi = pclosestApproach.GetPhi(eyeCS)/degree;
              if (thePhi > -1 && thePhi < 181) {

                const double depth = depthProfile.Y(distance);

                if (xUpGeom < 0 || depth > xUpGeom)
                  xUpGeom = depth;
                if (xLowGeom < 0 || depth < xLowGeom)
                  xLowGeom = depth;

                if (xUpDAQ < 0 || depth > xUpDAQ)
                  xUpDAQ = depth;
                if (xLowDAQ < 0 || depth < xLowDAQ)
                  xLowDAQ = depth;

              }
            }
          }
        }
      } else {
        if (fVerbosity > 10)
          cout << " skipping mirror with less than "
               << minSLTPixels
               << " pixels " << candidatePixelIds.size() << endl;
      }
    }
  } catch (OutOfBoundException& e) {
    // ignore
  }

  if (xUpGeom > 0 && xLowGeom > 0) {
    xUpGeom /= (g/cm2);
    xLowGeom /= (g/cm2);
  }

  if (xUpDAQ > 0 && xLowDAQ > 0) {
    xUpDAQ /= (g/cm2);
    xLowDAQ /= (g/cm2);
  }

  if (fVerbosity > -1) {
    cout << "geom FOV boundaries: " << detEye.GetId() << " "
         << xUpGeom << " " << xLowGeom << "\n"
            "DAQ  FOV boundaries: " << detEye.GetId() << " "
         << xUpDAQ << " " << xLowDAQ << endl;
    if (xUpGeom != xUpDAQ || xLowGeom != xLowDAQ)
      cout << " ######### GEOM-FOV update! ##########" << endl;
    if (fabs(xUp - xUpDAQ) > 10*g/cm2 || fabs(xLow - xLowDAQ) > 10*g/cm2)
      cout << " ######### FOV update! ##########" << endl;
  }

  xLow = xLowDAQ;
  xUp = xUpDAQ;

  return true;
}


bool
FOVCalculator::FdUpAndRunning()
  const
{
  bool oneIsUp = false;
  const fdet::FDetector& fdDet = det::Detector::GetInstance().GetFDetector();
  try {
    for (auto eyeIter = fdDet.EyesBegin(fdet::FDetComponentSelector::eAll);
         eyeIter != fdDet.EyesEnd(fdet::FDetComponentSelector::eAll); ++eyeIter) {
      for (auto telIter = eyeIter->TelescopesBegin();
           telIter != eyeIter->TelescopesEnd(); ++telIter) {
        if (telIter->GetDAQStatus()) {
          oneIsUp = true;
          break;
        }
      }
      if (oneIsUp)
        break;
    }
  } catch (NonExistentComponentException& c) {
    cerr << " no FdUptime available!!\n"
         << c.GetExceptionName() << "\n " << c.GetMessage() << endl;
  }

  return oneIsUp;
}