FdApertureLightFinderKG/FdApertureLightFinder.cc

Go to the documentation of this file.

#include "FdApertureLightFinder.h"

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

#include <utl/Reader.h>
#include <utl/ErrorLogger.h>
#include <utl/MathConstants.h>
#include <utl/Vector.h>
#include <utl/PhysicalConstants.h>
#include <utl/TabulatedFunctionErrors.h>

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

#include <fevt/FEvent.h>
#include <fevt/Eye.h>
#include <fevt/Telescope.h>
#include <fevt/EyeHeader.h>
#include <fevt/EyeRecData.h>
#include <fevt/TelescopeRecData.h>
#include <fevt/PixelRecData.h>
#include <fevt/Pixel.h>
#include <fevt/FdConstants.h>

#include <det/Detector.h>

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

#include <utl/Vector.h>

#include <vector>
#include <list>
#include <algorithm>
#include <iostream>
#include <sstream>

using namespace FdApertureLightFinderKG;
using namespace fwk;
using namespace utl;
using namespace evt;
using namespace fevt;
using namespace det;
using namespace std;


fwk::VModule::ResultFlag
FdApertureLightFinder::Init()
{
  CentralConfig* const cc = CentralConfig::GetInstance();

  Branch topB = cc->GetTopBranch("FdApertureLightKG");

  if (!topB) {
    ERROR("Could not find branch FdApertureLightKG");
    return eFailure;
  }

  topB.GetChild("verbosity").GetData(fVerbosity);
  topB.GetChild("onlyUsedForTimeFit").GetData(fOnlyUsedForTimeFit);

  Branch zetaB = topB.GetChild("zetaOptimization");
  // min and max zeta angle
  zetaB.GetChild("minZetaAngle").GetData(fminZetaAngle);
  zetaB.GetChild("maxZetaAngle").GetData(fmaxZetaAngle);

  // zeta search step width
  zetaB.GetChild("stepZetaAngle").GetData(fstepZetaAngle);

  // safety margin to enlarge S/N max result
  zetaB.GetChild("safetyMargin").GetData(fSafetyMargin);

  return eSuccess;
}


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

  int nGoodEyes = 0;
  for (FEvent::EyeIterator eyeiter = event.GetFEvent().EyesBegin(ComponentSelector::eHasData);
       eyeiter != event.GetFEvent().EyesEnd(ComponentSelector::eHasData); ++eyeiter) {

    Eye& eye = *eyeiter;

    if (CalculateShowerGeometryData(eye)) {
      if (FindLightFlux(eye)) {
        if (fVerbosity >= 0) {
          ostringstream msg;
          msg << "Calculated telescope- and eye-based light at aperture for eye " << eye.GetId();
          INFO(msg);
        }
        ++nGoodEyes;
      } else {
        if (fVerbosity >= 0) {
          ostringstream msg;
          msg << "Failed to calculate light at aperture for eye " << eye.GetId();
          INFO(msg);
        }
      }
    }

  }

  if (!nGoodEyes) {
    if (fVerbosity > 0)
      INFO("No good eyes, returning eContinueLoop");
    return eContinueLoop;
  }

  return eSuccess;
}


bool
FdApertureLightFinder::CalculateShowerGeometryData(const Eye& eye)
{
  fShowerGeometryData.clear();

  if (!eye.HasRecData() || !eye.GetRecData().HasFRecShower()) {
    ostringstream msg;
    msg << "No RecData for eye " << eye.GetId() << " skipping...";
    if (fVerbosity >= 0)
      INFO(msg);
    return false;
  }

  const fdet::FDetector& detFD   = det::Detector::GetInstance().GetFDetector();
  const fdet::Eye& detEye        = detFD.GetEye(eye);
  CoordinateSystemPtr eyeCS      = detEye.GetEyeCoordinateSystem();

  //---> get reconstructed axis and time fit parameters

  const ShowerFRecData& shower   = eye.GetRecData().GetFRecShower();

  const utl::Vector& axis        = shower.GetAxis(); // upward!
  const utl::Point& core         = shower.GetCorePosition();
  const utl::TimeStamp& coreTime = shower.GetCoreTime();

  const TimeStamp& eyeGPSTime    = eye.GetHeader().GetTimeStamp();
  const TimeStamp eyeTriggerTime = eyeGPSTime-TimeInterval(1000*100*utl::ns);

  //---> get time of first and last pulsed pixel

  double minT_i  = 1e9;
  double maxT_i  = 0;
  bool firstTime = true;

  for (fevt::EyeRecData::ConstPixelIterator pixiter = eye.GetRecData().TimeFitPixelsBegin();
       pixiter != eye.GetRecData().TimeFitPixelsEnd();
       ++pixiter) {

    if (not pixiter->HasRecData())
        continue;

    const fdet::Channel& detChannel = detFD.GetChannel(*pixiter);
    const double timebinsize = detChannel.GetFADCBinSize();

    const unsigned int timeoffset = eye.GetTelescope( pixiter->GetTelescopeId() ).GetTimeOffset();
    const double t1 = timebinsize * pixiter->GetRecData().GetPulseStart() + timeoffset;
    const double t2 = timebinsize * pixiter->GetRecData().GetPulseStop() + timeoffset;

    if (t2 > maxT_i || firstTime)
      maxT_i = t2;
    if (t1 < minT_i || firstTime)
      minT_i = t1;

    if (firstTime)
      firstTime = false;
  }

  if (maxT_i <= minT_i) {
    INFO("shower time length is <= 0, skipping...");
    return false;
  }

  //---> fill chi_i vector and corresponding ADC bin for each telescope

  for (fevt::Eye::ConstTelescopeIterator teliter = eye.TelescopesBegin(ComponentSelector::eHasData);
       teliter != eye.TelescopesEnd(ComponentSelector::eHasData); ++teliter) {
    const fdet::Telescope& detTel = detFD.GetTelescope(*teliter);

    const unsigned int telId       = teliter->GetId();
    const fdet::Camera& detCam     = detTel.GetCamera();
    const double timeBinSize       = detCam.GetFADCBinSize();
    const unsigned int adcTraceLength = detCam.GetFADCTraceLength();
    const unsigned int timeoffset  = teliter->GetTimeOffset();
    const utl::Vector& telAxis     = detTel.GetAxis();

    // for distance from shower core to emission point (aDist)
    // law of cosines: bDist^2 = aDist^2 + cDist^2 - 2*aDist*cDist*cos(beta)
    // cDist = distance(telescope, core)
    // bDist-aDist = (telTime-coreTime)/c = delta (see below)
    const utl::Point& telPos       = detTel.GetPosition();
    const Vector coreTelescopeVec  = telPos-core;
    const double cDist             = coreTelescopeVec.GetMag();
    const double cosBeta           = CosAngle(axis, coreTelescopeVec);

    // The localTelCS is at the tel. pos., but unlike the TelescopeCoordinateSystem,
    // its Z-axis points towards the zenith and not towards the aperture.
    // This is intended as an equivalent of the eyeCS, but at the telescope.
    // Unfortunately, a telescope has no backwall angle, so we're using the site's
    // X/Y directions.
    CoordinateSystemPtr localTelCS = fwk::LocalCoordinateSystem::Create(telPos);

    vector<ShowerGeometryData> telData;

//#warning This is a stop-gap measure to prevent aperture-light being calculated on bins that were used for the baseline
    const unsigned int nBaselineBins = 200;// FIXME constant should come from FdCalibrator
    for (unsigned int i = nBaselineBins; i < adcTraceLength; ++i) {

      const double t1    =   i   * timeBinSize + timeoffset;
      const double t2    = (i+1) * timeBinSize + timeoffset;
      const double tMean = t1 + 0.5 * timeBinSize;

      // "mean" below refering to the middle of the time bin as with tMean
      utl::TimeStamp meanTimeAtTelescope = eyeTriggerTime + tMean;

      const double deltaMean =
        (meanTimeAtTelescope - coreTime).GetInterval()*kSpeedOfLight;
      const double meanDistance =
        -(cDist*cDist - deltaMean*deltaMean) / (2*(deltaMean + cDist*cosBeta));
      const utl::Point meanPointOnShower = core - meanDistance*axis;

      // Vector from telescope to point on shower axis that corresponds to chi_i.
      utl::Vector chi_i_vec = meanPointOnShower - telPos;
      chi_i_vec.Normalize();

      // outside FOV wrt. the telescopes with margin for zeta
      // Note: This could be refined by checking chi_i_vec's vertical and horizontal
      //       angles to the telescope axis separately because the telescope is rectagular
      //       and not circular.
      const double angleToTelAxis = Angle(chi_i_vec, telAxis);
      if (angleToTelAxis > 25*deg)
        continue;

      ShowerGeometryData tmpShGeoData;
      tmpShGeoData.fADCTraceIndex    = i;
      tmpShGeoData.fT1               = t1;
      tmpShGeoData.fT2               = t2;
      tmpShGeoData.fChiVecMean[0]    = chi_i_vec.GetX(localTelCS);
      tmpShGeoData.fChiVecMean[1]    = chi_i_vec.GetY(localTelCS);
      tmpShGeoData.fChiVecMean[2]    = chi_i_vec.GetZ(localTelCS);
      tmpShGeoData.fTelFluxSum       = 0;
      tmpShGeoData.fTelVarianceSum   = 0;
      tmpShGeoData.fTelBgVarianceSum = 0;
      tmpShGeoData.fEyeFluxSum       = 0;
      tmpShGeoData.fEyeVarianceSum   = 0;
      tmpShGeoData.fEyeBgVarianceSum = 0;
      tmpShGeoData.fNearBorder       = false;
      tmpShGeoData.fUsedForTimeFit   = (t1 >= minT_i && t2 <= maxT_i);

      telData.push_back(tmpShGeoData);
    }

    fShowerGeometryData[telId] = telData;
  }

  return true;
}


bool
FdApertureLightFinder::FindLightFlux(Eye& eye)
{
  // rough stepping in 0.5 degree steps
  double zStep = 0.5 * degree;
  double zMin  = fminZetaAngle;
  double zMax  = fmaxZetaAngle + zStep;
  double zeta  = FindZeta(eye, zMin, zMax, zStep);

  // restepping around max
  zMin  = zeta - zStep;
  zMax  = zeta + zStep;
  zStep = fstepZetaAngle;
  if (zMin < 0.1*degree)
    zMin = 0.1*degree;

  zeta = FindZeta(eye, zMin, zMax, zStep);

  if (zeta < 0) {
    INFO("Zeta search failed - Bailing out for this eye");
    return false;
  }

  // enlarge zeta to collect more light
  zeta += fSafetyMargin;

  // Since we're fitting a eye-global zeta for now,
  // we save the same zeta in all telescopes and the eye
  eye.GetRecData().SetZeta(zeta);
  for (fevt::Eye::TelescopeIterator zetaTelIter = eye.TelescopesBegin(ComponentSelector::eHasData);
       zetaTelIter != eye.TelescopesEnd(ComponentSelector::eHasData); ++zetaTelIter) {
    if (!zetaTelIter->HasRecData())
      zetaTelIter->MakeRecData();
    zetaTelIter->GetRecData().SetZeta(zeta);
  }

  // find and fill telescope light fluxes at maximum
  FindZeta(eye, zeta, zeta, zStep, true);

  // order telescopes in time
  vector<LastTelTime> timeOrderedTelescopes;
  for (fevt::Eye::ConstTelescopeIterator teliter = eye.TelescopesBegin(ComponentSelector::eHasData);
       teliter != eye.TelescopesEnd(ComponentSelector::eHasData); ++teliter) {
    const unsigned int telId = teliter->GetId();
    LastTelTime telTime(telId, 0); // just a struct of two unsigned ints for sorting
    timeOrderedTelescopes.push_back(telTime);

    const vector<ShowerGeometryData>& telescopeData = fShowerGeometryData[telId];

    // backwards iteration for finding the last time with signal
    for (int i = telescopeData.size()-1; i >= 0; --i) {
      if (telescopeData[i].fEyeVarianceSum > 0) {
        timeOrderedTelescopes.back().fLastTime = telescopeData[i].fT2;
        break;
      }
    }
  } // end for each telescope
  sort(timeOrderedTelescopes.begin(), timeOrderedTelescopes.end(), LastTelTime::TelescopeTimeCmp);

  vector<double> telFluoTime;
  vector<double> telFluoTimeBin;
  vector<double> telFluoFlux;
  vector<double> telVFluoFlux;
  vector<double> telVBackground;
  vector<const vector<unsigned int>*> telPixelsInZeta; // Pointers on purpose. STL Optimization!
  vector<double> eyeFluoTime;
  vector<double> eyeFluoTimeBin;
  vector<double> eyeFluoFlux;
  vector<double> eyeVFluoFlux;
  vector<double> eyeVBackground;

  // loop over telescopes ordered in time
  unsigned int latestFluxTime = 0;
  for (vector<LastTelTime>::const_iterator telescopeIter = timeOrderedTelescopes.begin();
       telescopeIter != timeOrderedTelescopes.end(); ++telescopeIter) {
    const int telId = telescopeIter->fTelId;

    const vector<ShowerGeometryData>& telescopeData = fShowerGeometryData[telId];
    typedef pair<double, double> FarFromBorderTimeRange;
    list<FarFromBorderTimeRange> spotFarFromBorderRanges;

    bool isNearBorder = true;
    double rangeStartTime = -1e99;

    for (unsigned int i = 0; i < telescopeData.size(); ++i) {
      const ShowerGeometryData& telDataNow = telescopeData[i];
      const double dT = telDataNow.fT2 - telDataNow.fT1;
      const double t_i = telDataNow.fT1 + 0.5*dT;

      if (isNearBorder) {
        if (!telDataNow.fNearBorder && telDataNow.fTelVarianceSum > 0) {
          // start a new "far from border" time range
          isNearBorder = false;
          rangeStartTime = telDataNow.fT1; // start of time bin, not center
        }
      } else {
        if (telDataNow.fNearBorder || !(telDataNow.fTelVarianceSum > 0) ) {
          // stop the current "far from border" time range
          isNearBorder = true;
          // runs until end of bin (T2), not center.
          // now we are near border => end of previous bin (T2) = start of current bin (T1)
          spotFarFromBorderRanges.push_back(std::make_pair(rangeStartTime, telDataNow.fT1));
          if (fVerbosity > 0) {
            ostringstream msg;
            msg << "Telescope " << telId << ": Detected new time range for "
                << "light-collection-efficiency-less profile reconstruction: "
                << rangeStartTime << " ns to " << telDataNow.fT1 << " ns.";
            INFO(msg);
          }
          rangeStartTime = -1e99;
        }
      }

      if (telDataNow.fTelVarianceSum > 0) {
        telFluoTime.push_back(t_i);
        telFluoTimeBin.push_back(dT);
        telFluoFlux.push_back(telDataNow.fTelFluxSum);
        telVFluoFlux.push_back(telDataNow.fTelVarianceSum);
        telVBackground.push_back(telDataNow.fTelBgVarianceSum);
        telPixelsInZeta.push_back(&(telDataNow.fPixelsInZeta));
      }

      if (telDataNow.fT1 > latestFluxTime && telDataNow.fUsedForTimeFit) {
        if (telDataNow.fEyeVarianceSum > 0) {
          eyeFluoTime.push_back(t_i);
          eyeFluoTimeBin.push_back(dT);
          eyeFluoFlux.push_back(telDataNow.fEyeFluxSum);
          eyeVFluoFlux.push_back(telDataNow.fEyeVarianceSum);
          eyeVBackground.push_back(telDataNow.fEyeBgVarianceSum);
        }
      }

    } // end for each time bin

    // If (for some obscure reason), we're still within a "far from border" range at
    // the end of the track, add another range to the list. This isn't known to happen. Ever.
    if (rangeStartTime >= 0 && !isNearBorder) {
      const ShowerGeometryData& last = telescopeData.back();
      spotFarFromBorderRanges.push_back(std::make_pair(rangeStartTime, last.fT2));
      if (fVerbosity > 0) {
        ostringstream msg;
        msg << "Telescope " << telId << ": Detected new time range for "
               "light-collection-efficiency-less profile reconstruction: "
            << rangeStartTime << " ns to " << last.fT2 << " ns. "
               " This is not expected for normal FD events. BE CAREFUL.";
        WARNING(msg);
      }
    }

    if (!telFluoTime.empty()) {
      Telescope& telescope = eye.GetTelescope(telId);
      FillTelescopeFluxes(telescope, telFluoTime, telFluoTimeBin, telFluoFlux,
                          telVFluoFlux, telVBackground, telPixelsInZeta,
                          spotFarFromBorderRanges);
    }

    // clear the telescope profile containers for reuse for the next telescope
    telFluoTime.clear();
    telFluoTimeBin.clear();
    telFluoFlux.clear();
    telVFluoFlux.clear();
    telVBackground.clear();
    telPixelsInZeta.clear();
    latestFluxTime = telescopeIter->fLastTime;
  } // end for each telescope

  if (eyeFluoTime.empty())
    return false;

  FillEyeFluxes(eye, eyeFluoTime, eyeFluoTimeBin, eyeFluoFlux, eyeVFluoFlux, eyeVBackground);
  return true;
}


TabulatedFunctionErrors&
FdApertureLightFinder::GetClearedTelescopeLightFlux(Telescope& tel,
                                                    const FdConstants::LightSource source)
{
  if (!tel.GetRecData().HasLightFlux(source))
    tel.GetRecData().MakeLightFlux(source);

  TabulatedFunctionErrors& flux = tel.GetRecData().GetLightFlux(source);
  flux.Clear();
  return flux;
}


TabulatedFunctionErrors&
FdApertureLightFinder::GetClearedEyeLightFlux(Eye& eye,
                                              const FdConstants::LightSource source)
{
  if (!eye.GetRecData().HasLightFlux(source))
    eye.GetRecData().MakeLightFlux(source);

  TabulatedFunctionErrors& flux = eye.GetRecData().GetLightFlux(source);
  flux.Clear();
  return flux;
}


// Fetches the data structures from the eye and calls the generic
// profile filler
void
FdApertureLightFinder::FillEyeFluxes(Eye& eye,
                                     const vector<double>& fluoTime,
                                     const vector<double>& fluoTimeBin,
                                     const vector<double>& fluoFlux,
                                     const vector<double>& vfluoFlux,
                                     const vector<double>& vBackground)
{
  TabulatedFunctionErrors& fluolightprofile =
    GetClearedEyeLightFlux(eye, FdConstants::eFluorTotal);
  TabulatedFunctionErrors& lightprofile =
    GetClearedEyeLightFlux(eye, FdConstants::eTotal);
  TabulatedFunctionErrors& backgroundProfile =
    GetClearedEyeLightFlux(eye, FdConstants::eBackground);

  FillFluxes(fluolightprofile, lightprofile, backgroundProfile,
             fluoTime, fluoTimeBin, fluoFlux, vfluoFlux, vBackground);
}


// Fetches the data structures from the telescope and calls the generic
// profile filler
void
FdApertureLightFinder::FillTelescopeFluxes(Telescope& tel,
                                           const vector<double>& fluoTime,
                                           const vector<double>& fluoTimeBin,
                                           const vector<double>& fluoFlux,
                                           const vector<double>& vfluoFlux,
                                           const vector<double>& vBackground,
                                           const vector<const vector<unsigned int>*>& pixelsInZeta,
                                           const list<pair<double, double>>& farFromBorderTimeRanges)
{
  TabulatedFunctionErrors& fluolightprofile =
    GetClearedTelescopeLightFlux(tel, FdConstants::eFluorTotal);
  TabulatedFunctionErrors& lightprofile =
    GetClearedTelescopeLightFlux(tel, FdConstants::eTotal);
  TabulatedFunctionErrors& backgroundProfile =
    GetClearedTelescopeLightFlux(tel, FdConstants::eBackground);

  FillFluxes( fluolightprofile, lightprofile, backgroundProfile,
              fluoTime, fluoTimeBin, fluoFlux, vfluoFlux, vBackground );

  tel.GetRecData().SetSpotFarFromBorderTimeRanges(farFromBorderTimeRanges);

  // Fill the "pixels in zeta at given time" information as well
  const unsigned int nBins = pixelsInZeta.size();
  vector<vector<unsigned int>>& zetaPixels = tel.GetRecData().GetPixelsInZetaOverTime();
  zetaPixels.resize(nBins);
  for (unsigned int iTimeBin = 0; iTimeBin < nBins; ++iTimeBin)
    zetaPixels[iTimeBin] = *pixelsInZeta[iTimeBin];
}


// The generic profile filler guts, called only from the eye and tel specific fillers
void
FdApertureLightFinder::FillFluxes(TabulatedFunctionErrors& fluolightprofile,
                                  TabulatedFunctionErrors& lightprofile,
                                  TabulatedFunctionErrors& backgroundProfile,
                                  const vector<double>& fluoTime,
                                  const vector<double>& fluoTimeBin,
                                  const vector<double>& fluoFlux,
                                  const vector<double>& vfluoFlux,
                                  const vector<double>& vBackground)
{
  for (unsigned int i = 0, n = fluoTime.size(); i < n; ++i) {
    const double sigmaFluoFlux = sqrt(vfluoFlux[i]);
    const double tMean = fluoTime[i];
    const double dt = fluoTimeBin[i];
    lightprofile.PushBack(tMean, 0.5*dt, fluoFlux[i], sigmaFluoFlux);
    fluolightprofile.PushBack(tMean, 0.5*dt, fluoFlux[i], sigmaFluoFlux);
    backgroundProfile.PushBack(tMean, 0.5*dt, 0, sqrt(vBackground[i]));
  }
}


double
FdApertureLightFinder::FindZeta(fevt::Eye& eye,
                                double initZeta,
                                double finZeta,
                                double stepZeta,
                                bool fillFlux)
{
  double maxSN = 0;
  double bestZeta = initZeta;

  if (fillFlux && (initZeta + 1e-12 < finZeta || initZeta - 1e-12 > finZeta)) {
    ERROR("You called FindZeta with the fillFlux=true argument, but the initial and final guesses "
          "for zeta are not identical. This is an error. Either call FindZeta for zeta optimization "
          "(fillFlux=false) or for filling the light fluxes (fillFlux=true, initial zeta==final zeta)");
    return -9999*degree;
  }

  const fdet::FDetector& detFD = det::Detector::GetInstance().GetFDetector();
  //CoordinateSystemPtr eyeCS = detFD.GetEye(eye).GetEyeCoordinateSystem();
  const double referenceLambda = detFD.GetReferenceLambda();

  // loop on zeta

  for (double zeta = initZeta; zeta <= finZeta; zeta += stepZeta) {

    const double cosZeta = cos(zeta);
    double totalSignal = 0;
    double totalNoise2 = 0;

    for (fevt::Eye::TelescopeIterator teliter = eye.TelescopesBegin(ComponentSelector::eHasData);
         teliter != eye.TelescopesEnd(ComponentSelector::eHasData); ++teliter) {

      const unsigned int telId = teliter->GetId();

      const fdet::Telescope& dettel =
        det::Detector::GetInstance().GetFDetector().GetTelescope(*teliter);
      CoordinateSystemPtr localTelCS = fwk::LocalCoordinateSystem::Create(dettel.GetPosition());
      const double gainVariance = dettel.GetCamera().GetGainVariance();
      const double area = dettel.GetDiaphragmArea();
      const double area2 = area*area;

      vector<ShowerGeometryData>& telescopeData = fShowerGeometryData[telId];

      for (fevt::Telescope::ConstPixelIterator pixiter = teliter->PixelsBegin(ComponentSelector::eHasData);
           pixiter != teliter->PixelsEnd(ComponentSelector::eHasData); ++pixiter) {

        if (!pixiter->HasRecData())
          continue;

        const fdet::Pixel& detPixel = detFD.GetPixel(*pixiter);
        // get photon-to-photo electron-factor
        const double photon2Pe      = detPixel.GetDiaPhoton2PEFactor(referenceLambda);

        const Vector& pixeldir      = detPixel.GetDirection();
        const double pixelX         = pixeldir.GetX(localTelCS);
        const double pixelY         = pixeldir.GetY(localTelCS);
        const double pixelZ         = pixeldir.GetZ(localTelCS);

        // photon trace as written by calibrator: (raw trace - baseline)*calibconst
        const TraceD& photontrace = pixiter->GetRecData().GetPhotonTrace();
        const double RMS2 = pow(pixiter->GetRecData().GetRMS(), 2);

        for (unsigned int i = 0; i < telescopeData.size(); ++i) {
          ShowerGeometryData& thisTelData = telescopeData[i];
          // zeta-optimization done only with data that was used for the time fit:
          if (!thisTelData.fUsedForTimeFit && !fillFlux)
            continue;

          const double cosAngle =   thisTelData.fChiVecMean[0] * pixelX
                                  + thisTelData.fChiVecMean[1] * pixelY
                                  + thisTelData.fChiVecMean[2] * pixelZ;

          if (cosAngle > cosZeta) {
            const unsigned int idx = thisTelData.fADCTraceIndex;
            if (idx >= photontrace.GetStart() && idx <= photontrace.GetStop()) {
              const double photonSignal_i = photontrace[idx];
              const double noise2_i =
                photonSignal_i < 0 ? RMS2 : RMS2 + photonSignal_i / photon2Pe * (1 + gainVariance);

              totalSignal += photonSignal_i;
              totalNoise2 += noise2_i;

              if (fillFlux) {
                thisTelData.fTelFluxSum       += photonSignal_i/area;
                thisTelData.fTelVarianceSum   += noise2_i/area2;
                thisTelData.fTelBgVarianceSum += RMS2/area2;
                // Fill the vector of pixels that are within Zeta at the given time (bin i)
                thisTelData.fPixelsInZeta.push_back(pixiter->GetId());

                if (!IsNearBorder(Vector(thisTelData.fChiVecMean[0],
                                         thisTelData.fChiVecMean[1],
                                         thisTelData.fChiVecMean[2],
                                         localTelCS),
                                  *teliter, cosZeta)) {
                  // if NOT near border
                  if (thisTelData.fUsedForTimeFit) {
                    thisTelData.fEyeFluxSum       += photonSignal_i/area;
                    thisTelData.fEyeVarianceSum   += noise2_i/area2;
                    thisTelData.fEyeBgVarianceSum += RMS2/area2;
                  } else if (fOnlyUsedForTimeFit) {
                    thisTelData.fNearBorder = true;
                  }
                } else { // NEAR border
                  thisTelData.fNearBorder = true;
                }
              }
            } // end for each charge
          } // end if pixel is in zeta
        } // time-loop
      } // pixel-loop
    } // telescope-loop

    const double signalToNoise = totalNoise2 > 0 ? totalSignal / sqrt(totalNoise2) : 0.;

    if (!fillFlux)
      cout << "            "
           << " zeta [deg] " << zeta/degree
           << " S/N " << signalToNoise << " S " << totalSignal
           << endl;

    if (signalToNoise > maxSN) {
      maxSN = signalToNoise;
      bestZeta = zeta;
    }
  } // zeta-loop

  if (maxSN < 1)
    bestZeta = -9999*degree;

  if (!fillFlux)
    cout << "            Best Zeta " << bestZeta/degree << "\n" << endl;

  return bestZeta;
}


bool
FdApertureLightFinder::IsNearBorder(const utl::Vector& direction,
                                    const fevt::Telescope& telescope,
                                    const double cosZeta)
  const
{
  for (fevt::Telescope::ConstMirrorEventBorderPixelsIterator iter = telescope.MirrorEventBorderPixelsBegin();
       iter != telescope.MirrorEventBorderPixelsEnd(); ++iter) {
    if (CosAngle(direction, *iter) > cosZeta)
      return true;
  }

  return false;
}