FDetector/Telescope.cc

Go to the documentation of this file.

#include <boost/lexical_cast.hpp>

#include <utl/Point.h>
#include <utl/Vector.h>
#include <utl/ErrorLogger.h>
#include <utl/AugerUnits.h>
#include <utl/Reader.h>
#include <utl/TimeStamp.h>
#include <utl/TimeInterval.h>
#include <utl/UTCDateTime.h>
#include <utl/TabulatedFunction.h>
#include <utl/MathConstants.h>
#include <utl/Math.h>
#include <utl/AugerException.h>
#include <utl/AugerUnits.h>
#include <utl/Md5Sum.h>
#include <utl/UTMPoint.h>
#include <utl/ReferenceEllipsoid.h>

#include <det/VManager.h>
#include <det/Detector.h>

#include <fevt/Pixel.h>

#include <fdet/FDetector.h>
#include <fdet/Channel.h>
#include <fdet/Pixel.h>
#include <fdet/Eye.h>
#include <fdet/Diaphragm.h>
#include <fdet/Mirror.h>
#include <fdet/Filter.h>
#include <fdet/Corrector.h>
#include <fdet/Camera.h>
#include <fdet/Telescope.h>

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

#include <set>

using namespace std;
using namespace utl;
using namespace det;
using namespace fwk;


namespace fdet {

  Telescope::Telescope(const unsigned int eyeId, const unsigned int telescopeId) :
    fEyeId(eyeId),
    fTelescopeId(telescopeId)
  {
    fEyeIdString = boost::lexical_cast<string>(fEyeId);
    fTelescopeIdString = boost::lexical_cast<string>(fTelescopeId);
  }


  Telescope::~Telescope()
  {
    for (InternalChannelCollection::const_iterator cIt = fChannel.begin();
         cIt != fChannel.end(); ++cIt)
      delete cIt->second;

    for (InternalPixelCollection::const_iterator pIt = fPixel.begin();
         pIt != fPixel.end(); ++pIt)
      delete pIt->second;

    delete fConfigSignatureStr;
    delete fConfigSignature;
    delete fPosition;
    delete fAxis;
    delete fCommissionTime;
    delete fDecommissionTime;
    delete fMeasuredRelativeEfficiency;
    delete fTelPointingId;
    delete fTelPointingPhiMap;
    delete fTelPointingElevationMap;
    delete fDiaphragm;
    delete fMirror;
    delete fFilter;
    delete fCorrector;
    delete fCamera;
  }


  const std::string&
  Telescope::GetConfigSignature(const std::string& module)
    const
  {
    if (module != fConfigSignatureModule)
      GetConfigSignatureStr(module);
    if (!fConfigSignature)
      fConfigSignature = new string(Md5Sum(GetConfigSignatureStr(module)).GetHexDigest());
    return *fConfigSignature;
  }


  const std::string&
  Telescope::GetConfigSignatureStr(const std::string& module)
    const
  {
    if (module != fConfigSignatureModule) {
      delete fConfigSignature;
      fConfigSignature = 0;
      delete fConfigSignatureStr;
      fConfigSignatureStr = 0;
      fConfigSignatureModule = module;
    }

    if (fConfigSignatureStr)
      return *fConfigSignatureStr;

    ostringstream configSS;

    configSS << "Module       : " << module << "\n";
    const double refWl = Detector::GetInstance().GetFDetector().GetReferenceLambda();
    configSS << "RefWavelength: " << refWl/nanometer << "nm\n";

    configSS << "Filter       : T=";
    bool first = true;
    const TabulatedFunction& filtTrans = GetFilter().GetTransmittance();
    for (unsigned int iF = 0; iF < filtTrans.GetNPoints(); ++iF)
      if (filtTrans[iF].Y() != 0) {
        configSS << (!first ? ", " : "") << filtTrans[iF].X()/nanometer << "nm:" << filtTrans[iF].Y();
        first = false;
      }
    configSS << '\n';

    const bool hasCorrector = HasCorrectorRing();
    configSS << "Corrector    : " << (hasCorrector ? "yes" : "no");
    if (hasCorrector) {
      configSS << '\n'
               << "               R1=" << GetCorrector().GetInnerRadius()/m
               << "m,  R2=" << GetCorrector().GetOuterRadius()/m
               << "m, Tm=" << GetCorrector().GetMeanLensThickness()/m << "m,\n"
               << "               RI=";
      first = true;
      const TabulatedFunction& lensRI = GetCorrector().GetRefractiveIndex();
      for (unsigned int iR = 0; iR<lensRI.GetNPoints(); ++iR) {
        if (lensRI[iR].Y() != 0) {
          configSS << (!first ? ", " : "") << lensRI[iR].X()/nanometer << "nm:" << lensRI[iR].Y();
          first = false;
        }
      }
      configSS << ",\n                T=";
      first = true;
      const TabulatedFunction& lensT = GetCorrector().GetTransmittance();
      for (unsigned int iT = 0; iT<lensT.GetNPoints(); ++iT) {
        if (lensT[iT].Y() != 0) {
          configSS << (!first ? ", " : "") << lensT[iT].X()/nanometer << "nm:" << lensT[iT].Y();
          first = false;
        }
      }
    }
    configSS << '\n';

    configSS << "Mirror       : Radius=" << GetMirror().GetRadiusOfCurvature()/m << "m,\n";
    configSS << "               R=";
    first = true;
    const TabulatedFunction& mirRef = GetMirror().GetReflectivity();
    for (unsigned int iM = 0; iM<mirRef.GetNPoints(); ++iM)
      if (mirRef[iM].Y() != 0) {
        configSS << (!first ? ", " : "") << mirRef[iM].X()/nanometer << "nm:" << mirRef[iM].Y();
        first = false;
      }
    configSS << '\n';

    const fdet::Pixel& detPixel = GetPixel(1);
    const fdet::Channel& detChannel = GetChannel(detPixel.GetChannelId()); // mapped channel
    configSS << "Camera       : Merc=" << GetCamera().GetMercedesEfficiency()
             << ", R_f=" << GetCamera().GetRadiusFocal()
             << ", H_m=" << GetCamera().GetMercedesHeight()
             << ", dd_m=" << GetCamera().GetMercedesBase()
             << ", nCols=" << GetCamera().GetLastColumn()-GetCamera().GetFirstColumn() + 1
             << ", nRows=" << GetCamera().GetLastRow()-GetCamera().GetFirstRow() + 1
             << ", Omega=" << GetCamera().GetEta()
             << ", gain=" << detChannel.GetElectronicsGain()
             << '\n';
    const TabulatedFunction& qEff = detPixel.GetQEfficiency();
    configSS << "               Qeff=";
    first = true;
    for (unsigned int iF = 0; iF < qEff.GetNPoints(); ++iF) {
      if (qEff[iF].Y() != 0) {
        configSS << (!first ? ", " : "") << qEff[iF].X()/nanometer << "nm:" << qEff[iF].Y();
        first = false;
      }
    }
    fConfigSignatureStr = new string(configSS.str());
    return *fConfigSignatureStr;
  }


  const Channel&
  Telescope::GetChannel(const unsigned int channelId)
    const
  {
    InternalChannelCollection::const_iterator cIt = fChannel.find(channelId);

    if (cIt != fChannel.end())

      return *cIt->second;

    else {
      if (channelId < GetFirstChannelId() || channelId > GetLastChannelId()) {
        ostringstream err;
        err << "Channel Id=" << channelId << " is out of range; "
               "min=" << GetFirstChannelId() << " max=" << GetLastChannelId();
        ERROR(err);
        throw NonExistentComponentException(err.str());
      }

      const Channel* const ch = new Channel(fEyeId, fTelescopeId,
                                            channelId,
                                            GetParentPhysicalEyeIdString(),
                                            GetParentPhysicalIdString());
      fChannel[channelId] = ch;
      return *ch;
    }
  }


  const Channel&
  Telescope::GetChannel(const fdet::Pixel& p)
    const
  {
    return GetChannel(p.GetChannelId());
  }


  const Channel&
  Telescope::GetChannel(const fevt::Pixel& p)
    const
  {
    return GetChannel(GetPixel(p.GetId()).GetChannelId());
  }


  Point
  Telescope::GetPosition()
    const
  {
    if (!fPosition) {
      const fdet::FDetector& theFDet = Detector::GetInstance().GetFDetector();
      const fdet::Eye& realEye = theFDet.GetEye(GetParentPhysicalEyeId());
      const CoordinateSystemPtr realEyeCoordinateSystem = realEye.GetEyeCoordinateSystem();
      double bay_angle;
      GetTelescopeData(bay_angle, "bay_angle", "telescope", "telescope bay angle");
      double bay_distance;
      GetTelescopeData(bay_distance, "bay_distance", "telescope", "telescope bay distance");
      fPosition = new Point(bay_distance*cos(bay_angle), bay_distance*sin(bay_angle), 0, realEyeCoordinateSystem);
    }
    return *fPosition;
  }


  CoordinateSystemPtr
  Telescope::GetTelescopeCoordinateSystem()
    const
  {
    if (!fCoordinateSystem) {
      const CoordinateSystemPtr telCS = LocalCoordinateSystem::Create(GetPosition());

      const double axisphi = GetAxis().GetPhi(telCS);

      const CoordinateSystemPtr auxCS =
        CoordinateSystem::RotationZ(axisphi, telCS);

      const double theta = GetAxis().GetTheta(auxCS);

      // Set z-axis of the telescope CS = telescope axis
      fCoordinateSystem = CoordinateSystem::RotationY(theta, auxCS);
    }
    return fCoordinateSystem;
  }


  bool
  Telescope::HasCorrectorRing()
    const
  {
    return GetDiaphragm().HasCorrectorRing();
  }


  double
  Telescope::GetDiaphragmRadius()
    const
  {
    return GetDiaphragm().GetRadius();
  }


  double
  Telescope::GetDiaphragmArea()
    const
  {
    return GetDiaphragm().GetArea();
  }


  const TimeStamp&
  Telescope::GetCommissionTime()
    const
  {
    if (!fCommissionTime) {
      string timeString;
      GetTelescopeData(timeString, "commission", "telescope", "telescope commission time");
      fCommissionTime = new TimeStamp(boost::lexical_cast<UTCDateTime>(timeString).GetTimeStamp());
    }

    return *fCommissionTime;
  }


  const TimeStamp&
  Telescope::GetDecommissionTime()
    const
  {
    if (!fDecommissionTime) {
      string timeString;
      GetTelescopeData(timeString, "decommission", "telescope", "telescope decommission time");
      fDecommissionTime = new TimeStamp(boost::lexical_cast<UTCDateTime>(timeString).GetTimeStamp());
    }

    return *fDecommissionTime;
  }


  const Pixel&
  Telescope::GetPixel(const unsigned int pixelId)
    const
  {
    if (pixelId < GetFirstPixelId() || pixelId > GetLastPixelId()) {
      ostringstream err;
      err << "Pixel Id=" << pixelId << " is out of range; "
          << "min=" << GetFirstPixelId() << " max=" << GetLastPixelId();
      ERROR(err);
      throw NonExistentComponentException(err.str());
    }

    const InternalPixelCollection::const_iterator pIt = fPixel.find(pixelId);

    if (pIt != fPixel.end())
      return *pIt->second;
    else {
      Pixel* const pix = new Pixel(fEyeId, fTelescopeId, pixelId,
                                   GetParentPhysicalEyeIdString(),
                                   GetParentPhysicalIdString());
      fPixel[pixelId] = pix;
      return *pix;
    }
  }


  const Pixel&
  Telescope::GetPixel(const unsigned int row, const unsigned int col)
    const
  {
    const unsigned int pixelId = GetLastRow()*(col-1) + row;
    return GetPixel(pixelId);
  }


  const Diaphragm&
  Telescope::GetDiaphragm()
    const
  {
    if (!fDiaphragm)
      fDiaphragm = new Diaphragm(fEyeId, fTelescopeId,
                                 GetParentPhysicalEyeIdString(),
                                 GetParentPhysicalIdString());
    return *fDiaphragm;
  }


  const Mirror&
  Telescope::GetMirror()
    const
  {
    if (!fMirror)
      fMirror = new Mirror(fEyeId, fTelescopeId,
                           GetParentPhysicalEyeIdString(),
                           GetParentPhysicalIdString());
    return *fMirror;
  }


  const Filter&
  Telescope::GetFilter()
    const
  {
    if (!fFilter)
      fFilter = new Filter(fEyeId, fTelescopeId,
                           GetParentPhysicalEyeIdString(),
                           GetParentPhysicalIdString());
    return *fFilter;
  }


  const Corrector&
  Telescope::GetCorrector()
    const
  {
    if (!fCorrector)
      fCorrector = new Corrector(fEyeId, fTelescopeId,
                                 GetParentPhysicalEyeIdString(),
                                 GetParentPhysicalIdString());
    return *fCorrector;
  }


  const Camera&
  Telescope::GetCamera()
    const
  {
    if (!fCamera)
      fCamera = new Camera(fEyeId, fTelescopeId,
                           GetParentPhysicalEyeIdString(),
                           GetParentPhysicalIdString());
    return *fCamera;
  }


  void
  Telescope::CacheTelescopePointing()
    const
  {
    if (!fTelPointingId)
      GetTelescopeData(fTelPointingId, "pointing", "telescope", "telescope optical axis default pointing id");

    if (!fTelPointingPhiMap) {
      GetTelescopeData(fTelPointingPhiMap, "opticalAxisPhi", "telescope", "telescope optical axis phi");

      if (!fTelPointingElevationMap)
        GetTelescopeData(fTelPointingElevationMap, "opticalAxisElevation", "telescope", "telescope optical axis elevation");

      if (IsVirtual()) {
        // The axis phi and elevation are defined relative to the eye coordinate
        // system. Here the values are transformed from physical eye coordinate
        // system to the one of the virtual eye.
        const fdet::FDetector& theFDet = Detector::GetInstance().GetFDetector();
        const fdet::Eye& eye = theFDet.GetEye(fEyeId);
        const fdet::Eye& realEye = theFDet.GetEye(GetParentPhysicalEyeId());
        const CoordinateSystemPtr eyeCoordinateSystem = eye.GetEyeCoordinateSystem();
        const CoordinateSystemPtr realEyeCoordinateSystem = realEye.GetEyeCoordinateSystem();

        for (PointingAngleSet::iterator phiIt = fTelPointingPhiMap->begin(),
            end = fTelPointingPhiMap->end(); phiIt != end; ++phiIt) {
          const PointingAngleSet::iterator elevationIt = fTelPointingElevationMap->find(phiIt->first);
          if (elevationIt == fTelPointingElevationMap->end())
            cout << "This should never happen! Found azimuth without corresponding elevation" << endl;

          const double azimuth = phiIt->second;
          const double elevation = elevationIt->second;

          utl::Vector axis(1., kPi/2 - elevation, azimuth,
                           realEyeCoordinateSystem, Vector::kSpherical);
          const double newAzimuth = axis.GetPhi(eyeCoordinateSystem);
          const double newElevation = kPi/2 - axis.GetTheta(eyeCoordinateSystem);

          phiIt->second = newAzimuth;
          elevationIt->second = newElevation;
        }
      } // end if IsVirtual
    }
  }


  std::string
  Telescope::GetTelescopePointingId()
    const
  {
    if (!fTelPointingId)
      CacheTelescopePointing();
    return *fTelPointingId;
  }


  Telescope::PointingAngleSet
  Telescope::GetTelescopePointingPhi()
    const
  {
    if (!fTelPointingPhiMap)
      CacheTelescopePointing();
    return *fTelPointingPhiMap;
  }


  Telescope::PointingAngleSet
  Telescope::GetTelescopePointingElevation()
    const
  {
    if (!fTelPointingElevationMap)
      CacheTelescopePointing();
    return *fTelPointingElevationMap;
  }


  utl::Vector
  Telescope::GetAxis()
    const
  {
    if (!fTelPointingId) {

      CacheTelescopePointing();

      if (fTelPointingElevationMap->find(*fTelPointingId) == fTelPointingElevationMap->end() ||
          fTelPointingPhiMap->find(*fTelPointingId) == fTelPointingPhiMap->end()) {

        ostringstream err;
        err << "Failure reading Fd telescope pointing (eye: " << fEyeId << ", tel: " << fTelescopeId << ") from telescope list XML: "
               "the pointing that corresponds to the Id: \'" << *fTelPointingId << "\' does not exist. The config knows about: \n";
        for (PointingAngleSet::const_iterator iPointing = fTelPointingElevationMap->begin();
             iPointing != fTelPointingElevationMap->end(); ++iPointing)
          err << " pointing id: \'" << iPointing->first << "\' at elevation: " << iPointing->second << "\n";
        ERROR(err);
        throw NonExistentComponentException(err.str());
      }
    } // end if we haven't read the pointing info yet

    if (!fAxis) {
      const double azimuth = fTelPointingPhiMap->find(*fTelPointingId)->second;
      const double elevation = fTelPointingElevationMap->find(*fTelPointingId)->second;

      const fdet::Eye& eye = Detector::GetInstance().GetFDetector().GetEye(fEyeId);
      const CoordinateSystemPtr eyeCoordinateSystem = eye.GetEyeCoordinateSystem();
      fAxis = new Vector(1, kPi/2 - elevation, azimuth,
                         eyeCoordinateSystem, Vector::kSpherical);
    }

    return *fAxis;
  }


  void
  Telescope::GenerateOutOfBorderPixels()
    const
  {
    // top and bottom
    const int rowInc = GetLastRow() - GetFirstRow() + 2;
    for (unsigned int row = GetFirstRow() - 1; row <= GetLastRow() + 1; row += rowInc) {
      // starting from first-1 and going to last+1 to include the corners
      for (unsigned int col = GetFirstColumn()-1; col <= GetLastColumn()+1; ++col) {
        const utl::Vector dir = GetCameraPixelDirection(row, col);
        fOutBorderPixelsDir.push_back(dir);
      }
    }

    // left and right
    const int colInc = GetLastColumn() - GetFirstColumn() + 2;
    for (unsigned int col = GetFirstColumn() - 1; col <= GetLastColumn() + 1;
         col += colInc) {
      for (unsigned int row = GetFirstRow(); row <= GetLastRow(); ++row) {
        const utl::Vector dir = GetCameraPixelDirection(row, col);
        fOutBorderPixelsDir.push_back(dir);
      }
    }
  }


  Telescope::OutOfBorderPixelsIterator
  Telescope::OutOfBorderPixelsBegin()
    const
  {
    if (!fOutBorderPixelsDir.size())
      GenerateOutOfBorderPixels();
    return fOutBorderPixelsDir.begin();
  }


  Telescope::OutOfBorderPixelsIterator
  Telescope::OutOfBorderPixelsEnd()
    const
  {
    if (!fOutBorderPixelsDir.size())
      GenerateOutOfBorderPixels();
    return fOutBorderPixelsDir.end();
  }


  OFFLINE_DEFINE_CONST_ITERATOR_RANGE(OutOfBorderPixelsIterator, Telescope, OutOfBorderPixels)


  double
  Telescope::GetModelWavelengthDependence(const double wavelength)
    const
  {
    const utl::TabulatedFunction& mirrorRefl = GetMirror().GetReflectivity();
    const utl::TabulatedFunction& filterTrans = GetFilter().GetTransmittance();

    // mercedes reflections, when it potentially gets wavelength-dependent ...
    const double mercEff = GetCamera().GetMercedesEfficiency();
    const double meanRefl = GetCamera().GetMercedesReflections();
    const double meanMercEff = pow(mercEff, meanRefl);

    const TabulatedFunction& Qeff = GetPixel(1).GetQEfficiency();

    double eff = mirrorRefl.Y(wavelength) * filterTrans.Y(wavelength) *
      Qeff.Y(wavelength) *
      meanMercEff; // for eventual future use

    if (HasCorrectorRing()) {
      const utl::TabulatedFunction& lensTrans = GetCorrector().GetTransmittance();
      const utl::TabulatedFunction& lensRI = GetCorrector().GetRefractiveIndex();
      const double meanLensThickness = GetCorrector().GetMeanLensThickness();

      //# internal transmission thru 11mm, x=transmission thru 10mm

      if (lensTrans.Y(wavelength) == 0)
        return 0;

      const double internal_absorption
        = exp(meanLensThickness/(10.*mm)*log(lensTrans.Y(wavelength)));

      //# transmission at air - glass interface, normal incidence
      // there is about 3.5% reflection from 0 to >20 deg
      const double n = lensRI.Y(wavelength);
      const double T_air_glass = 1 - Sqr((n-1)/(n+1));

      //# fraction of aperture covered by corrector ring
      // effectively fracArea is larger, since shadowing by the camera is mostly in the hole of the lens
      const double fracArea = 1 - Sqr(0.85/1.1);

      //# total effective transmittance (plus 2 glass-air interfaces)
      const double lensTransmission = 1 + fracArea * (internal_absorption * Sqr(T_air_glass) - 1);
      //# just internal transmittance
      //const double lensTransmission = 1 + fracArea * (internal_absorption - 1);

      eff *= lensTransmission;
    }
    return eff;
  }


  double
  Telescope::GetModelRelativeEfficiency(const double wavelength)
    const
  {
    const double refWl = Detector::GetInstance().GetFDetector().GetReferenceLambda();
    return GetModelWavelengthDependence(wavelength) / GetModelWavelengthDependence(refWl);
  }


  double
  Telescope::GetModelMeanEfficiency(const string& configSignature, const double wavelength)
    const
  {
    const CoordinateSystemPtr& telCS = GetTelescopeCoordinateSystem();
    double meanPixelEff = 0;
    int nPix = 0;
    for (unsigned int i = GetFirstPixelId(); i <= GetLastPixelId(); ++i) {
      ++nPix;
      const Pixel& pixel = GetPixel(i);
      const double cosTheta = pixel.GetDirection().GetCosTheta(telCS);
      meanPixelEff += 1 / (pixel.GetSimulatedEndToEndCalibration(configSignature) * cosTheta);
    }
    meanPixelEff *= GetModelRelativeEfficiency(wavelength);
    if (nPix)
      meanPixelEff /= nPix;
    return meanPixelEff;
  }


  double
  Telescope::GetMeasuredRelativeEfficiency(const double wavelength)
    const
  {
    const utl::TabulatedFunction& eff = GetMeasuredRelativeEfficiency();
    return eff.Y(wavelength);
  }


  const utl::TabulatedFunction&
  Telescope::GetMeasuredRelativeEfficiency()
    const
  {
    if (!fMeasuredRelativeEfficiency) {
      utl::TabulatedFunction measuredEff;
      GetTelescopeData(measuredEff, "relativeResponse", "telescope", "telescope relative response");
      const double refWl = Detector::GetInstance().GetFDetector().GetReferenceLambda();
      fMeasuredRelativeEfficiency = new TabulatedFunction;
      for (unsigned int i = 0; i < measuredEff.GetNPoints(); ++i) {
        fMeasuredRelativeEfficiency->PushBack(measuredEff.GetX(i), measuredEff.GetY(i)/measuredEff.Y(refWl));
      }
    }
    return *fMeasuredRelativeEfficiency;
  }


  // delegate these issues to the camera
  unsigned int
  Telescope::GetFirstRow()
    const
  {
    return GetCamera().GetFirstRow();
  }


  unsigned int
  Telescope::GetFirstColumn()
    const
  {
    return GetCamera().GetFirstColumn();
  }


  unsigned int
  Telescope::GetLastRow()
    const
  {
    return GetCamera().GetLastRow();
  }


  unsigned int
  Telescope::GetLastColumn()
    const
  {
    return GetCamera().GetLastColumn();
  }


  unsigned int
  Telescope::GetLastPixelId()
    const
  {
    const Camera& cam = GetCamera();
    return cam.GetLastColumn() * cam.GetLastRow();
  }


  double
  Telescope::GetUpTimeFraction()
    const
  {
    UpdateFdUpTime();
    return fUpTimeFraction;
  }


  bool
  Telescope::IsInAquisition()
    const
  {
    UpdateFdUpTime();
    return fUpTimeFraction > 0;
  }


  int
  Telescope::GetDAQStatus()
    const
  {
    UpdateFdUpTime();
    return fStatus;
  }


  void
  Telescope::UpdateFdUpTime()
    const
  {
    if (fUpTimeValidityStamp.IsValid())
      return;

    int startTime;
    GetTelescopeData(startTime, "uptime_gpsStart", "telescope", "start validity telescope up time");
    int endTime;
    GetTelescopeData(endTime, "uptime_gpsStop", "telescope", "end validity telescope up time");
    GetTelescopeData(fUpTimeFraction, "uptime_fraction", "telescope", "telescope up time fraction");
    GetTelescopeData(fStatus, "status", "telescope", "telescope DAQ status");
    fUpTimeValidityStamp.SetValidityInterval(startTime, endTime);
  }


  void
  Telescope::Update()
  {
    // Update any volatile telescope properties

    if (fIdOfParentPhysicalEye < 0) {
      fIdOfParentPhysicalEye = 0; // default to being a physical telescope
      fIdOfParentPhysicalTelescope = 0;
      GetTelescopeData(fIdOfParentPhysicalEye, "fromEye", "virtualTelescope",
                       "parent eye id of virtual telescope");
      GetTelescopeData(fIdOfParentPhysicalTelescope, "fromTelescope", "virtualTelescope",
                       "parent telescope id of virtual telescope");

      fIdOfParentPhysicalEyeString
        = boost::lexical_cast<string>(fIdOfParentPhysicalEye);
      fIdOfParentPhysicalTelescopeString
        = boost::lexical_cast<string>(fIdOfParentPhysicalTelescope);
    }

    delete fTelPointingId;
    fTelPointingId = 0;

    delete fAxis;
    fAxis = 0;

    delete fConfigSignature;
    fConfigSignature = 0;
    delete fConfigSignatureStr;
    fConfigSignatureStr = 0;

    fCoordinateSystem.reset();
    fOutBorderPixelsDir.clear();

    if (fDiaphragm)
      fDiaphragm->Update();
    if (fMirror)
      fMirror->Update();
    if (fFilter)
      fFilter->Update();
    if (fCorrector)
      fCorrector->Update();
    if (fCamera)
      fCamera->Update();

    for (InternalChannelCollection::const_iterator cIt = fChannel.begin();
         cIt != fChannel.end(); ++cIt)
      cIt->second->Update();

    for (InternalPixelCollection::const_iterator pIt = fPixel.begin();
         pIt != fPixel.end(); ++pIt)
      pIt->second->Update();
  }


  void
  Telescope::CachePixelCalibrations()
    const
  {
    if (fCalibValidityStamp.IsValid())
      return;

    vector<int> pixelStatus;
    GetTelescopeData(pixelStatus, "status", "fd_calib", "pixel status");
    const unsigned int statusSize = pixelStatus.size();

    // Delete all pixel calibs in the telescope if
    // the validty interval has expired and
    // get new calibration status
    for (unsigned int i = GetFirstPixelId(); i <= GetLastPixelId(); ++i) {
      const Pixel& pixel = GetPixel(i);
      pixel.fWavelengthCalib.Clear();
      if (i >= statusSize)
        pixel.fPixelStatus = Pixel::eUnknown;
      else
        switch (pixelStatus[i]) {
        case 0:
          pixel.fPixelStatus = Pixel::eGood;
          break;
        case 1:
          pixel.fPixelStatus = Pixel::eBadCalibration;
          break;
        case 2:
          pixel.fPixelStatus = Pixel::eUnknown;
          break;
        default:
          pixel.fPixelStatus = Pixel::eUnknown;
          break;
        }
    }

    GetTelescopeData(fCalibMap, "pixel_calibs", "fd_calib", "pixel calibration");

    // Fill pixels from data retrieved above
    for (unsigned int iPixelId = GetFirstPixelId();
         iPixelId <= GetLastPixelId(); ++iPixelId) {

      if (fCalibMap.find(iPixelId) != fCalibMap.end()) {

        TabulatedFunction& pixCal = GetPixel(iPixelId).fWavelengthCalib;

        // done in this tortuous way because TabulatedFunction is so crap and
        // has no sensible assignment operator
        TabulatedFunction& cachedTabFunc = fCalibMap.find(iPixelId)->second;

        for (TabulatedFunction::Iterator tIt = cachedTabFunc.Begin();
             tIt != cachedTabFunc.End(); ++tIt) {
          pixCal.PushBack(tIt->X(), tIt->Y());
        }

        pixCal.SetBoundaryTolerance(0.1*nanometer);

      } else {

        const string err = "Failure reading FD calibration from database";
        ERROR(err);
        throw NonExistentComponentException(err);

      }
    }

    int startTime;
    GetTelescopeData(startTime, "start_time", "fd_calib", "pixel calibration start validity");
    int endTime;
    GetTelescopeData(endTime, "end_time", "fd_calib", "pixel calibration end validity");
    fCalibValidityStamp.SetValidityInterval(startTime, endTime);
  }

  void
  Telescope::CachePixelOpticalEfficiencyCorrections()
    const
  {
    if (fOpticalEfficiencyValidityStamp.IsValid())
      return;

    vector<int> pixelStatus;
    GetTelescopeData(pixelStatus, "status", "fd_optical_efficiency", "pixel status");
    const unsigned int statusSize = pixelStatus.size();

    // Delete all pixel optical efficiency corrections in the telescope if
    // the validty interval has expired and
    // get new optical efficiency correction status
    for (unsigned int i = GetFirstPixelId(); i <= GetLastPixelId(); ++i) {
      const Pixel& pixel = GetPixel(i);
      pixel.fWavelengthOpticalEfficiencyCorrection.Clear();
      if (i >= statusSize)
        pixel.fPixelOpticalEfficiencyStatus = Pixel::eUnknown;
      else
        switch (pixelStatus[i]) {
        case 0:
          pixel.fPixelOpticalEfficiencyStatus = Pixel::eGood;
          break;
        case 1:
          pixel.fPixelOpticalEfficiencyStatus = Pixel::eBadCalibration;
          break;
        case 2:
          pixel.fPixelOpticalEfficiencyStatus = Pixel::eUnknown;
          break;
        default:
          pixel.fPixelOpticalEfficiencyStatus = Pixel::eUnknown;
          break;
        }
    }

    GetTelescopeData(fOpticalEfficiencyMap, "pixel_optical_efficiency_corrections", "fd_optical_efficiency", "pixel optical efficiency correction");

    // Fill pixels from data retrieved above
    for (unsigned int iPixelId = GetFirstPixelId();
         iPixelId <= GetLastPixelId(); ++iPixelId) {

      if (fOpticalEfficiencyMap.find(iPixelId) != fOpticalEfficiencyMap.end()) {

        TabulatedFunction& pixOpticalEfficiency = GetPixel(iPixelId).fWavelengthOpticalEfficiencyCorrection;

        // done in this tortuous way because TabulatedFunction is so crap and
        // has no sensible assignment operator
        TabulatedFunction& cachedTabFunc = fOpticalEfficiencyMap.find(iPixelId)->second;

        for (TabulatedFunction::Iterator tIt = cachedTabFunc.Begin();
             tIt != cachedTabFunc.End(); ++tIt) {
          pixOpticalEfficiency.PushBack(tIt->X(), tIt->Y());
        }

        pixOpticalEfficiency.SetBoundaryTolerance(0.1*nanometer);

      } else {

        const string err = "Failure reading FD optical efficiency correction from database";
        ERROR(err);
        throw NonExistentComponentException(err);

      }
    }

    int startTime;
    GetTelescopeData(startTime, "start_time", "fd_optical_efficiency", "pixel optical efficiency start validity");
    int endTime;
    GetTelescopeData(endTime, "end_time", "fd_optical_efficiency", "pixel optical efficiency end validity");
    fOpticalEfficiencyValidityStamp.SetValidityInterval(startTime, endTime);
  }

  bool
  Telescope::CachePixelCloudData()
    const
  {
    // Get all relevant telescope data from the cloud DB in a single shot.
    // The intention is that this function is not called during the normal
    // detector Update, but only when Pixel cloud data are explicitly requested

    try {
      const VManager& manager = Detector::GetInstance().
        GetAManagerRegister().GetManager("ACloudSQLManager");

      vector<vector<int> > pixel_coverage_start_end;

      const VManager::IndexMap indexMap({
        { "eyeId", GetParentPhysicalEyeIdString() },
        { "telescopeId", GetParentPhysicalIdString() }
      });

      const VManager::Status status =
        manager.GetData(pixel_coverage_start_end, "cloud_pixel", "pixel_id, cloud_coverage, start_time, end_time", indexMap);

      // Exception if no data are found.
      if (status == VManager::eNotFound || pixel_coverage_start_end.empty())
        return false;

      // check if all pixels are present
      set<int> pixelsInCloudDB;
      for (unsigned int i = 0; i < pixel_coverage_start_end.size(); ++i)
        pixelsInCloudDB.insert(pixel_coverage_start_end[i][0]);

      if ( pixelsInCloudDB.size() != GetLastPixelId()-GetFirstPixelId()+1 ) {
        ostringstream warn;
        warn << "no full telescope in cloud DB --> skip"
             << "(nPix=" << fPixel.size() << " nCloud="
             << pixelsInCloudDB.size() << ")";
        for (unsigned int i = 0; i < pixel_coverage_start_end.size(); ++i) {
          cout << " pixel Id: " << pixel_coverage_start_end[i][0]
               << " gpsStart: " << pixel_coverage_start_end[i][2]
               << " gpsStop: " << pixel_coverage_start_end[i][3]
               << endl;
        }
        WARNING(warn);
        return false;
      }

      // Unpack the cloud data from the manager.  Warning: No further error checks!
      for (unsigned int i = 0; i < pixel_coverage_start_end.size(); ++i) {
        const int pixelId = pixel_coverage_start_end[i][0];
        const int cloudIndex = pixel_coverage_start_end[i][1];
        const int gpsStart = pixel_coverage_start_end[i][2];
        const int gpsEnd = pixel_coverage_start_end[i][3];

        const Pixel& pixel = GetPixel(pixelId);
        pixel.fCloudIndex = cloudIndex;
        pixel.fCloudValidityStamp.SetValidityInterval(gpsStart, gpsEnd);
      }
    } catch (NonExistentComponentException&) {
      return false;
    }

    return true;
  }


  utl::Vector
  Telescope::GetCameraPixelDirection(const unsigned int pixelId)
    const
  {
    if (pixelId < GetFirstPixelId() ||
        pixelId > GetLastPixelId()) {
      ostringstream err;
      err << "Pixel number " << pixelId << " "
             "is not in the Camera";
      ERROR(err);
      throw NonExistentComponentException(err.str());
    }
    const fdet::Pixel& pix = GetPixel(pixelId);
    return GetCameraPixelDirection(pix.GetRow(), pix.GetColumn());
  }


  utl::Vector
  Telescope::GetCameraPixelDirection(const unsigned int row, const unsigned int col)
    const
  {
    // algorithm to generate the pixel directions in the camera system
    // Developer note:
    // If you ever change this, make sure it works for a row/column of pixels
    // *beyond* the camera as well. This is used for the out-of-border-pixels!
    const double eta0 = GetCamera().GetEta();
    const double centerRow = GetCamera().GetCenterRow();
    const double centerColumn = GetCamera().GetCenterColumn();
    const double dOmega = eta0;
    const double dPhi = kSqrt3 * eta0/2;

    // const double oo = (col - ((row%2) ? 10. : 10.5)) * dOmega;
    // const double ph = (11.66667 - row) * dPhi;

    const double oo = (col - centerColumn + (row%2)/2.) * dOmega;
    const double ph = (centerRow - row) * dPhi;

    const double mcosOO = -cos(oo);
    const double z = mcosOO * cos(ph);
    const double x = mcosOO * sin(ph);
    const double y = sin(oo);

    return Vector(-x, -y, -z, GetTelescopeCoordinateSystem());
  }


  double
  Telescope::GetModelMinWavelength()
    const
  {
    double minWavelength = GetMirror().GetReflectivity().XFront();
    minWavelength = max(minWavelength, GetFilter().GetTransmittance().XFront());
    minWavelength = max(minWavelength, GetCorrector().GetTransmittance().XFront());
    minWavelength = max(minWavelength, GetPixel(1).GetQEfficiency().XFront()); // alternatively: loop over all pixels
    GetMeasuredRelativeEfficiency(); // just to initialize fMeasuredEfficiency
    minWavelength = max(minWavelength, fMeasuredRelativeEfficiency->XFront());
    return minWavelength;
  }


  double
  Telescope::GetModelMaxWavelength()
    const
  {
    double maxWavelength = GetMirror().GetReflectivity().XBack();
    maxWavelength = min(maxWavelength, GetFilter().GetTransmittance().XBack());
    maxWavelength = min(maxWavelength, GetCorrector().GetTransmittance().XBack());
    maxWavelength = min(maxWavelength, GetPixel(1).GetQEfficiency().XBack()); // alternatively: loop over all pixels
    GetMeasuredRelativeEfficiency(); // just to initialize fMeasuredEfficiency
    maxWavelength = min(maxWavelength, fMeasuredRelativeEfficiency->XBack());
    return maxWavelength;
  }


  // ###########################################################################

  template<typename T>
  inline
  const T&
  Telescope::GetTelescopeData(T*& requestedData,
                              const string& property,
                              const string& component,
                              const string& errorMsg)
    const
  {
    if (!requestedData) {
      requestedData = new T;
      GetTelescopeData(*requestedData, property, component, errorMsg);
    }
    return *requestedData;
  }


  template<typename T>
  inline
  void
  Telescope::GetTelescopeData(T& requestedData,
                              const string& property,
                              const string& component,
                              const string& errorMsg)
    const
  {
    const VManager& manager = Detector::GetInstance().GetFManagerRegister();

    VManager::IndexMap indexMap;
    // If this is a virtual telescope , fetch the information from the
    // real telescope that this virtual telescope resembles
    if (IsVirtual() && component != "virtualTelescope") {
      indexMap["eyeId"] = fIdOfParentPhysicalEyeString;
      indexMap["telescopeId"] = fIdOfParentPhysicalTelescopeString;
    } else {
      indexMap["eyeId"] = fEyeIdString;
      indexMap["telescopeId"] = fTelescopeIdString;
    }

    const VManager::Status status =
      manager.GetData(requestedData, property, component, indexMap);

    if (status == VManager::eNotFound) {
      ostringstream err;
      err << "Did not find requested component for " << errorMsg << "; "
          << VManager::QueryInfoMessage(property, component, indexMap);
      ERROR(err);
      throw NonExistentComponentException(err.str());
    }
  }

}