FdApertureLightFinderOG/FdApertureLightFinder.cc

Go to the documentation of this file.

#include "FdApertureLightFinder.h"

#include <fwk/CentralConfig.h>
#include <fwk/CoordinateSystemRegistry.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/EyeRecData.h>
#include <fevt/PixelRecData.h>
#include <fevt/Pixel.h>
#include <fevt/FdComponentSelector.h>

#include <det/Detector.h>

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

#include <utl/Vector.h>

#include <sstream>
#include <iomanip>

using namespace FdApertureLightFinderOG;
using namespace fwk;
using namespace utl;
using namespace evt;
using namespace fevt;
using namespace det;
using namespace fdet;
using namespace std;

const double FdApertureLightFinder::kTelescopeMaxTheta = 89*utl::degree;
const double FdApertureLightFinder::kTelescopeMinTheta = 59*utl::degree;


fwk::VModule::ResultFlag
FdApertureLightFinder::Init(){

  Branch topB =
    CentralConfig::GetInstance()->GetTopBranch("FdApertureLightFinder");

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

  topB.GetChild("minZetaAngle").GetData(fminZetaAngle);
  topB.GetChild("maxZetaAngle").GetData(fmaxZetaAngle);
  topB.GetChild("stepZetaAngle").GetData(fstepZetaAngle);
  topB.GetChild("safetyMargin").GetData(fSafetyMargin);
  topB.GetChild("nEquivMin").GetData(fnEquivMin);
  topB.GetChild("nProfMax").GetData(fnProfMax);

  fTStart=0;
  fTStop =0;

  return eSuccess;
}


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

  FEvent& fevent = event.GetFEvent();

  FEvent::EyeIterator eyeiter;

  int nGoodEyes = 0;
  for (eyeiter = fevent.EyesBegin(ComponentSelector::eHasData);
       eyeiter != fevent.EyesEnd(ComponentSelector::eHasData);
       ++eyeiter){

    fevt::Eye& eye = *eyeiter;

    if (!eye.HasRecData()) {
      INFO("T0 must be taken from EyeRecData, skip this Eye");
      continue;
    }

    if (!eye.GetRecData().HasFRecShower()){
      INFO("No Geometrical info for this Eye - skip this Eye");
      continue;
    }

    if ( FindLightFlux(eye) )
      nGoodEyes ++;

  } // end eyeiter-loop

  if ( nGoodEyes == 0 )
    return eContinueLoop;
  else
    return eSuccess;

}


fwk::VModule::ResultFlag
FdApertureLightFinder::Finish()
{
  return eSuccess;
}


bool
FdApertureLightFinder::FindLightFlux(fevt::Eye& eye)
{

  // find first and last hit pixels to find shower start and stop
  if ( !FindSignalTimeRange(eye) )
    return false;

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

  eye.GetRecData().SetZeta(zeta);

  const fdet::FDetector& detFD = det::Detector::GetInstance().GetFDetector();
  const CoordinateSystemPtr& eyeCS = detFD.GetEye(eye).GetEyeCoordinateSystem();
  const Point& eyePos = detFD.GetEye(eye).GetPosition();
  const EyeRecData& eyerecdata =  eye.GetRecData();
  const double T0 = eyerecdata.GetTZero();
  const double chi0 = eyerecdata.GetChiZero();
  const double Rp = eyerecdata.GetRp();
  const Vector axis = eyerecdata.GetFRecShower().GetAxis();
  const Point core = eyerecdata.GetFRecShower().GetCorePosition();
  const Vector eyeToCore = core - eyePos;

//#warning fixed time binning at aperture !!!
  const double timeBinSize = 100*ns; // light at aperture binning fixed !!

  deque <double> fluoTime, fluoFlux, vfluoFlux, vBackground;

  // loop on time slots
  for (double time= fTStart;
       time <= fTStop;
       time += timeBinSize) {

    double flux_i=0;   //total light flux for that time
    double fluxVariance_i=0;
    double backgroundVariance_i=0;

    // vector from the detector to the light spot at start of that time bin
    const double chi_1 = chi0 - 2*atan( kSpeedOfLight/Rp *(time - T0) );
    double l = eyeToCore.GetMag() * sin(chi_1) / sin(chi0 - chi_1);
    if(l<0)
      l*=-1;
    Vector chi_i_vec1 = eyeToCore + l * axis;
    chi_i_vec1.Normalize();

    // vector from the detector to the light spot at end of that time bin
    const double chi_2 = chi0 - 2*atan( kSpeedOfLight/Rp *(time - T0+timeBinSize) );
    l = eyeToCore.GetMag() * sin(chi_2) / sin(chi0 - chi_2);
    if(l<0)
      l*=-1;
    Vector chi_i_vec2 = eyeToCore + l * axis;
    chi_i_vec2.Normalize();

    // vector from the detector to the light spot at mid of that time bin
    Vector chi_i_vec = chi_i_vec1 + chi_i_vec2;
    chi_i_vec.Normalize();

    if (chi_i_vec.GetTheta(eyeCS) > kTelescopeMaxTheta  ||
        chi_i_vec.GetTheta(eyeCS) < kTelescopeMinTheta)
      continue;

    unsigned int ntelescopes=0; // keep track if signal is from
                                // one or two telescopes

    const double referenceLambda = detFD.GetReferenceLambda();

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

      double telescope_flux_i=0;   // light flux seen by this telescope
      double telescope_variance_i=0;
      double telescope_bgVariance_i=0;

      double telescopeTimeOffset = teliter->GetTimeOffset();

      const fdet::Telescope& detTel = detFD.GetTelescope(*teliter);
      const double area        = detTel.GetDiaphragmArea();
      const double squaredArea = area*area;
      const double gainVariance = detTel.GetCamera().GetGainVariance();

      bool hasSignal = false;    // signal present?
      bool isNearBorder = false; // spot near the border?

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

        const fdet::Pixel& detPixel = detFD.GetPixel(*pixiter);
        const Vector& pixeldir = detPixel.GetDirection();
        const double photon2Pe = detPixel.GetDiaPhoton2PEFactor(referenceLambda);
        
        const double alpha=(Angle(chi_i_vec1, pixeldir) + Angle(chi_i_vec2, pixeldir))/2.;

        // check the pixel is inside the spot
        if(alpha < zeta) {
          
          hasSignal = true;
          isNearBorder = IsNearBorder(chi_i_vec, detTel,zeta);
          // if the spot is near the border, discard this time slot
          // for this telescope
          if (isNearBorder) {
            telescope_flux_i=0;
            telescope_variance_i=0;
            telescope_bgVariance_i=0;
            break;
          }

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

          const TraceD& photontrace= pixiter->GetRecData().GetPhotonTrace();

          int idx = int((time - telescopeTimeOffset)/timeBinSize);

          if(idx>=(int) photontrace.GetStart()&&
             idx<=(int) photontrace.GetStop()) {

            // cosangle correction not included
            // should be already in the end to end calibration
            telescope_flux_i += photontrace[idx]/area;

            const double baseLineVariance =
              pow(pixiter->GetRecData().GetRMS(),2)/squaredArea;

            telescope_bgVariance_i += baseLineVariance;

            if (photontrace[idx]<0)
              telescope_variance_i  += baseLineVariance;
            else {
              const double signalVariance =
                photontrace[idx]/photon2Pe * (1.+gainVariance);
              telescope_variance_i  +=
                (signalVariance/squaredArea + baseLineVariance);
            }
          }
        } // if angle
      }// for pixel


      if (!isNearBorder && hasSignal ){
        // only if the spot is totally contained add total flux
        flux_i+= telescope_flux_i;
        fluxVariance_i+= telescope_variance_i;
        backgroundVariance_i+= telescope_bgVariance_i;
        ntelescopes++;

      }  // if isfullycontained
    } // for telescope

    if (ntelescopes) {
      flux_i/=ntelescopes; // average in case signal is from two telescopes
      fluxVariance_i/=ntelescopes;
      backgroundVariance_i/=ntelescopes;
    }

    // fill deques with fluxes for subsequent combination
    if (flux_i){
      fluoTime.push_back(time);
      fluoFlux.push_back(flux_i);
      vfluoFlux.push_back(fluxVariance_i);
      vBackground.push_back(backgroundVariance_i);
    }

  } // for time


  // FdPulseFinder artifically enlarges pulse width
  // --> pop a few bins at the beginning and end
  //     of light flux until above 2*RMS
  unsigned int nPopFront=0;
  for ( unsigned int i=0;i<fluoTime.size();i++) {
    if(fluoFlux[i]>2.*sqrt(vfluoFlux[i]))
      break;
    nPopFront++;
  }
  for ( unsigned int i=0;i<nPopFront;i++) {
    fluoTime.pop_front();
    fluoFlux.pop_front();
    vfluoFlux.pop_front();
    vBackground.pop_front();
  }

  unsigned int nPopBack=0;
  for ( int i=(int) fluoTime.size()-1;i>-1;i--) {
    if(fluoFlux[i]>2.*sqrt(vfluoFlux[i]))
      break;
    nPopBack++;
  }
  for ( unsigned int i=0;i<nPopBack;i++) {
    fluoTime.pop_back();
    fluoFlux.pop_back();
    vfluoFlux.pop_back();
    vBackground.pop_back();
  }

  // combine bins and fill them in event structure
  if(fluoTime.size()<1)
    return false;
  else {
    CombineAndFillFluxes(eye, fluoTime, fluoFlux, vfluoFlux, vBackground);
    return true;
  }

}


//************************************************************************
bool
FdApertureLightFinder::FindSignalTimeRange(const fevt::Eye& eye)
{

  const fdet::FDetector& detFD = det::Detector::GetInstance().GetFDetector();
  const EyeRecData& eyerecdata =  eye.GetRecData();

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

  EyeRecData::ConstPixelIterator pixiter;
  for (pixiter = eyerecdata.TimeFitPixelsBegin();
       pixiter != eyerecdata.TimeFitPixelsEnd(); ++pixiter){

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

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

    const double timeoffset =
      eye.GetTelescope(pixiter->GetTelescopeId()).GetTimeOffset();
    const fevt::PixelRecData & pixelRecData = pixiter->GetRecData();
    const double t1 = timeBinSize*pixelRecData.GetPulseStart()+timeoffset;
    const double t2 = timeBinSize*pixelRecData.GetPulseStop()+timeoffset;

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

  }
  fTStart = minT_i;
  fTStop  = maxT_i;

  if (fTStop <= fTStart)
    return false;
  else
    return true;

}


//************************************************************************
void
FdApertureLightFinder::CombineAndFillFluxes(fevt::Eye& eye,
                                            const deque<double>& fluoTime,
                                            const deque<double>& fluoFlux,
                                            const deque<double>& vfluoFlux,
                                            const deque<double>& vBackgroundFlux)
{

//#warning Fixed 100ns time binning at the aperture
  const double timeBinSize = 100. * ns; // light at aperture binning fixed !!

  // make light fluxes in the event if necessary and fill
  // both total fluorescence light and total light.
  // During Cherenkov subtraction, fluorescence light will
  // undergo corrections

  if ( !eye.GetRecData().HasLightFlux(fevt::FdConstants::eFluorTotal) )
    eye.GetRecData().MakeLightFlux(fevt::FdConstants::eFluorTotal);

  if (!eye.GetRecData().HasLightFlux(fevt::FdConstants::eTotal) )
    eye.GetRecData().MakeLightFlux(fevt::FdConstants::eTotal);

  if (!eye.GetRecData().HasLightFlux(fevt::FdConstants::eBackground) )
    eye.GetRecData().MakeLightFlux(fevt::FdConstants::eBackground);

  TabulatedFunctionErrors& fluoLightProfile =
    eye.GetRecData().GetLightFlux(fevt::FdConstants::eFluorTotal);

  TabulatedFunctionErrors& lightProfile =
    eye.GetRecData().GetLightFlux(fevt::FdConstants::eTotal);

  TabulatedFunctionErrors& backgroundProfile =
    eye.GetRecData().GetLightFlux(fevt::FdConstants::eBackground);

  fluoLightProfile.Clear();
  lightProfile.Clear();
  backgroundProfile.Clear();

  // ------------  combine bins with low statistics --------------

  vector<double> c_fluoTime, c_fluoDTime, c_fluoFlux, c_vfluoFlux, c_vbgFlux;
  double nBins=0;
  double t1=fluoTime[0];
  double sum=0.;
  double sum2=0.;
  double bgSum2=0;
  double lastTime=fluoTime[0];

  //#warning here
  //const double timeBinSize = fdet::Channel::GetFADCBinSize();

  for ( unsigned int i=0;i<fluoTime.size();i++) {

    if(t1==0.)
      t1=fluoTime[i];

    // stop averaging if encounter a hole
    if(fluoTime[i]-lastTime>timeBinSize) {
      if(nBins>0.) {
        double dt    = fluoTime[i-1]-t1+timeBinSize;
        c_fluoTime.push_back(t1);
        c_fluoDTime.push_back(dt);
        c_fluoFlux.push_back(sum);
        c_vfluoFlux.push_back(sum2);
        c_vbgFlux.push_back(bgSum2);
        nBins=0.; sum=0.; sum2=0.; bgSum2=0.; t1=0.;
      }
    }

    if(t1==0.)
      t1=fluoTime[i];

    nBins++;
    sum += fluoFlux[i];
    sum2+= vfluoFlux[i];
    bgSum2+= vBackgroundFlux[i];
    lastTime=fluoTime[i];
    double nEquiv=sum*sum/sum2;
    if( fnEquivMin<=0 || nEquiv>fnEquivMin || i==fluoTime.size()-1 ) {
      double dt    = fluoTime[i]-t1+timeBinSize;
      c_fluoTime.push_back(t1);
      c_fluoDTime.push_back(dt);
      c_fluoFlux.push_back(sum);
      c_vfluoFlux.push_back(sum2);
      c_vbgFlux.push_back(bgSum2);
      nBins=0.; sum=0.; sum2=0.; bgSum2=0.;  t1=0.;
    }
  }

  // ------------  combine bins if above fnProfMax --------------
  if(c_fluoTime.size()>(unsigned int) fnProfMax) {

    double nCombi=(double) c_fluoTime.size()/ (double) fnProfMax;
    double combined=0.;

    t1=c_fluoTime[0];
    sum=0.;
    sum2=0.;
    bgSum2=0;
    lastTime=t1;

    for (unsigned int i=0;i<c_fluoTime.size();i++) {

      if(t1==0.)
        t1=c_fluoTime[i];

      // stop averaging if encounter a hole
      if(sum>0.&&c_fluoTime[i]-lastTime>c_fluoDTime[i-1]) {
        double dt    = c_fluoTime[i-1]-t1+c_fluoDTime[i-1];
        lightProfile.PushBack     (t1+dt/2.,dt/2.,sum,sqrt(sum2));
        fluoLightProfile.PushBack (t1+dt/2.,dt/2.,sum,sqrt(sum2));
        backgroundProfile.PushBack(t1+dt/2.,dt/2., 0.,sqrt(bgSum2));
        sum=0.; sum2=0.; t1=0.; bgSum2=0.; combined=0;
      }

      combined+=1.;
      sum += c_fluoFlux[i];
      sum2+= c_vfluoFlux[i];
      bgSum2+= c_vbgFlux[i];
      lastTime=c_fluoTime[i];

      if(t1==0.)
        t1=c_fluoTime[i];
      if(combined>=nCombi||i==c_fluoTime.size()-1) {
        double dt = c_fluoTime[i]-t1+c_fluoDTime[i];
        lightProfile.PushBack     (t1+dt/2.,dt/2.,sum,sqrt(sum2));
        fluoLightProfile.PushBack (t1+dt/2.,dt/2.,sum,sqrt(sum2));
        backgroundProfile.PushBack(t1+dt/2.,dt/2., 0.,sqrt(bgSum2));
        sum=0.; sum2=0.; t1=0.; bgSum2=0.; combined=0;
      }
    }
  }
  else {
    for (unsigned int i=0;i<c_fluoTime.size();i++) {
      sum    = c_fluoFlux[i];
      sum2   = c_vfluoFlux[i];
      bgSum2 = c_vbgFlux[i];
      double dt = c_fluoDTime[i];
      lightProfile.PushBack     (c_fluoTime[i]+dt/2.,dt/2.,sum,sqrt(sum2));
      fluoLightProfile.PushBack (c_fluoTime[i]+dt/2.,dt/2.,sum,sqrt(sum2));
      backgroundProfile.PushBack(c_fluoTime[i]+dt/2.,dt/2., 0.,sqrt(bgSum2));
    }
  }
}


bool
FdApertureLightFinder::IsNearBorder(const utl::Vector& direction,
                                    const fdet::Telescope& telescope,
                                    double zeta)
  const
{

  for (fdet::Telescope::OutOfBorderPixelsIterator iter =
         telescope.OutOfBorderPixelsBegin();
       iter!= telescope.OutOfBorderPixelsEnd();
       ++iter){

    if (Angle(direction,*iter) < zeta)
      return true;
  }

  return false;

}

unsigned int
FdApertureLightFinder::GetMaxOffset(const fevt::Eye& eye)
  const
{

  unsigned int maxoffset = 0;

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

    const unsigned int tmp = teliter->GetTimeOffset();

    if (tmp > maxoffset)
      maxoffset = tmp;

  }

  return maxoffset;
}


// New (faster) version of FindZeta by M. Unger
double
FdApertureLightFinder::FindZeta(const fevt::Eye& eye, double initZeta,
                                double finZeta, double stepZeta)
  const
{

  ostringstream info;
  info << "Eye: " << eye.GetId() << '\n';

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

  const EyeRecData& eyerecdata =  eye.GetRecData();
  const double T0 = eyerecdata.GetTZero();
  const double chi0 = eyerecdata.GetChiZero();
  const double Rp = eyerecdata.GetRp();
  const Vector axis = eyerecdata.GetFRecShower().GetAxis();
  const Point core = eyerecdata.GetFRecShower().GetCorePosition();
  const Vector eyeToCore = core - eyePos;

//#warning fixed time binning at aperture !!!
  const double timeBinSize = 100.*ns;

  double maxSN = 0;
  double bestZeta = initZeta;

  // M.U. : precalculate components of chi_i array for nested loop

  const int nTimes=(int) ((fTStop-fTStart)/timeBinSize)+1;
  // create vectors to hold component
  double * chi_i_vecX = new double[nTimes];
  double * chi_i_vecY = new double[nTimes];
  double * chi_i_vecZ = new double[nTimes];

  const double distEyeToCore = eyeToCore.GetMag();
  double time;
  int iTime=0;
  for (time= fTStart;
       time <= fTStop;
       time+= timeBinSize) {  // loop over time bins

    const double chi_i = chi0-2.* atan( kSpeedOfLight/Rp * (time - T0) );
    double l =  distEyeToCore * sin(chi_i) / sin(chi0 - chi_i);
    if(l<0)
      l*=-1;
    Vector chi_i_vec = eyeToCore + l * axis;
    chi_i_vec.Normalize();

    if (chi_i_vec.GetTheta(eyeCS) > kTelescopeMaxTheta  ||
        chi_i_vec.GetTheta(eyeCS) < kTelescopeMinTheta) {

      chi_i_vecX[iTime] = 9999;   // flag that outside telescope
      chi_i_vecY[iTime] = 9999;   // flag that outside telescope
      chi_i_vecZ[iTime] = 9999;   // flag that outside telescope

    } else {

      chi_i_vecX[iTime] = chi_i_vec.GetX(eyeCS);
      chi_i_vecY[iTime] = chi_i_vec.GetY(eyeCS);
      chi_i_vecZ[iTime] = chi_i_vec.GetZ(eyeCS);

    }

    iTime++;
  }

  // loop on zeta

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

    double signalOverNoise=0;
    double signalSum =0;
    double varianceSum =0;

    // M.U. : optimize speed by changing loop order

    double  charge_i=0.,noise2_i=0;

    const double referenceLambda = detFD.GetReferenceLambda();

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

      const unsigned int telescopeTimeOffset = teliter->GetTimeOffset();

      const fdet::Telescope& detTel = detFD.GetTelescope(*teliter);
      const double gainVariance = detTel.GetCamera().GetGainVariance();

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

        const fdet::Pixel& detPixel = detFD.GetPixel(*pixiter);
        const Vector& pixeldir = detPixel.GetDirection();
        const double photon2Pe = detPixel.GetDiaPhoton2PEFactor(referenceLambda);

        const double pixelX = pixeldir.GetX(eyeCS);
        const double pixelY = pixeldir.GetY(eyeCS);
        const double pixelZ = pixeldir.GetZ(eyeCS);

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

        const TraceD& photontrace = pixiter->GetRecData().GetPhotonTrace();
        const double baseLineVariance = pow(pixiter->GetRecData().GetRMS(),2);


        int iTime=0;
        for (time= fTStart;  // loop on time bin
             time <= fTStop;
             time+= timeBinSize) {

          if(*(chi_i_vecX+iTime)<100) { // inside telescope

            const double angle = acos(*(chi_i_vecX+iTime)*pixelX+
                                      *(chi_i_vecY+iTime)*pixelY+
                                      *(chi_i_vecZ+iTime)*pixelZ);


            if (angle<zeta) {
              const int idx = int((time - telescopeTimeOffset)/timeBinSize);

              if( idx>=(int) photontrace.GetStart()&&
                  idx<=(int) photontrace.GetStop()) {

                charge_i  = photontrace[idx];
                if ( charge_i < 0 )
                  noise2_i  = baseLineVariance;
                else {
                  const double signalVariance =
                    charge_i / photon2Pe * (1.+gainVariance);
                  noise2_i  = baseLineVariance + signalVariance;
                }

                signalSum += charge_i;
                varianceSum += noise2_i;
              }

            } //if zeta
          } // if telescope
          iTime++;

        }// for time
      } // for pixel
    } // for telescope


    if (varianceSum)
      signalOverNoise = signalSum / sqrt(varianceSum);

    info << "            "
         << " zeta [deg] " << setw(5) << zeta/degree
         << ", S/N " << setw(10) << signalOverNoise
         << ", S " << setw(9) << signalSum << '\n';

    if (signalOverNoise > maxSN) {
      maxSN = signalOverNoise;
      bestZeta = zeta;
    }
  } // for zeta

  delete []  chi_i_vecX;
  delete []  chi_i_vecY;
  delete []  chi_i_vecZ;

  if (maxSN < 1)
    bestZeta = -9999*degree;
  info << "            Best Zeta " << bestZeta/degree;
  INFO(info);

  return bestZeta;
}


// Configure (x)emacs for this file ...
// Local Variables:
// mode:c++
// compile-command: "make -C .. -k"
// End: