G/ApertureLight.cc

Go to the documentation of this file.

#include "ApertureLight.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 FdProfileConstrainedGeometryFitPG;
using namespace fwk;
using namespace utl;
using namespace evt;
using namespace fevt;
using namespace det;
using namespace std;


bool
ApertureLight::Init()
{
  Branch topB = fwk::CentralConfig::GetInstance()->GetTopBranch("FdProfileConstrainedGeometryFitPG");

  //topB.GetChild("verbosity").GetData(fVerbosity);
  fVerbosity = 0;

  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);

  // margin between shower position and non-recorded pixels
  zetaB.GetChild("borderMargin").GetData(fBorderMargin);

  fStripMode = false;

  return true;
}


bool
ApertureLight::Run(evt::Event& event)
{
  cout << "\n    --[PCGF internal ApertureLight]--> " << endl;

  if ( !event.HasFEvent() )
    return true;

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

  } // end for eyes

  if ( nGoodEyes == 0 ) {
    if (fVerbosity > 0)
      INFO("No good eyes, returning false");
    return false;
  }
  else
    return true;

}

bool
ApertureLight::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 int timeBinSize          = (int) detCam.GetFADCBinSize(); // FIXME, this should be a double!
    const 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; (int)i < adcTraceLength; ++i) {

      const int t1       =   i   * timeBinSize + timeoffset;
      const int t2       = (i+1) * timeBinSize + timeoffset;
      const double tMean = (double)t1 + 0.5 * (double)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 (!fStripMode && 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
ApertureLight::FindLightFlux(Eye& eye)
{
  // rough stepping in 0.5 degree steps
  double zStep = .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, telFluoTimeBin, telFluoFlux, telVFluoFlux, telVBackground;
  vector<const vector<unsigned int>*> telPixelsInZeta; // Pointers on purpose. STL Optimization!
  vector<double> eyeFluoTime, eyeFluoTimeBin, eyeFluoFlux, eyeVFluoFlux, 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 = -1.e99;

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

      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.
          //VN - mistake? - here we are NearBorder => far from border ends at fT1 of this bin
          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 = -1.e99;
        }
      }

      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.";
        INFO(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;
  else {
    FillEyeFluxes(eye, eyeFluoTime, eyeFluoTimeBin,
                  eyeFluoFlux, eyeVFluoFlux, eyeVBackground);
    return true;
  }
}


TabulatedFunctionErrors&
ApertureLight::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&
ApertureLight::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
ApertureLight::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
ApertureLight::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
ApertureLight::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; i < fluoTime.size(); ++i) {
    const double sigmaFluoFlux = sqrt(vfluoFlux[i]);
    const double tMean         = fluoTime[i];
    const double dt            = fluoTimeBin[i];
    lightprofile.PushBack     ( tMean, dt/2., fluoFlux[i], sigmaFluoFlux );
    fluolightprofile.PushBack ( tMean, dt/2., fluoFlux[i], sigmaFluoFlux );
    backgroundProfile.PushBack( tMean, dt/2., 0.,  sqrt(vBackground[i]) );
  }
}


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

  if (fillFlux && (initZeta+1.e-12 < finZeta || initZeta-1.e-12 > finZeta)) {
    ostringstream msg;
    msg << "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)";
    ERROR(msg);
    return -9999.*degree;
  }

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

  // fevt::EyeRecData& eyerecdata = eye.GetRecData();
  // const Vector& fSDP = eyerecdata.GetSDP();

  // 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 Vector& fSDP = teliter->GetRecData().GetSDP();

      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); // RMS2 = square of baseline RMS

        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 && (fStripMode || !fillFlux))
            continue;

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

          double angleSDP = Angle(fSDP, pixeldir);
//VN cout<<"SDP angle "<<angleSDP/deg<<endl;
          if (angleSDP > 90.*deg) {
            angleSDP = 180.*deg - angleSDP;
          }

          if ((!fStripMode && cosAngle > cosZeta) || (fStripMode && abs(angleSDP - 90.*deg) < zeta)) {
            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 (fStripMode || !IsNearBorder(Vector(thisTelData.fChiVecMean[0],
                                         thisTelData.fChiVecMean[1],
                                         thisTelData.fChiVecMean[2],
                                         localTelCS),
                                  *teliter, cosZeta))
                {
                  if (thisTelData.fUsedForTimeFit) {
                    thisTelData.fEyeFluxSum       += photonSignal_i/area;
                    thisTelData.fEyeVarianceSum   += noise2_i/area2;
                    thisTelData.fEyeBgVarianceSum += RMS2/area2;
                  }
                } // end if NOT near border
                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
ApertureLight::IsNearBorder(const utl::Vector& direction,
                                    const fevt::Telescope& telescope,
                                    double cosZeta)
  const
{
  double cosDeg = 0;
  if (fBorderMargin > 0) {
    cosDeg = cos(fBorderMargin);
  } else {
    cosDeg = cosZeta;
  }
  for (fevt::Telescope::ConstMirrorEventBorderPixelsIterator iter = telescope.MirrorEventBorderPixelsBegin();
       iter != telescope.MirrorEventBorderPixelsEnd(); ++iter) {
    if (CosAngle(direction, *iter) > cosDeg)
      return true;
  }
  return false;
}