UseMcGeometry.cc

Go to the documentation of this file.

#include "UseMcGeometry.h"

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

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

#include <evt/Event.h>

#include <sevt/SEvent.h>

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

#include <evt/ShowerSimData.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerFRecData.h>
#include <evt/ShowerSRecData.h>

#include <det/Detector.h>

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

#include <iostream>


using namespace std;
using namespace fwk;
using namespace utl;
using namespace fevt;
using namespace evt;
using namespace UseMcGeometryOG;



// --------------------------------------------------------------------------
VModule::ResultFlag
UseMcGeometry::Init()
{
  Branch topB =
    CentralConfig::GetInstance()->GetTopBranch("UseMcGeometry");

  topB.GetChild("setSDP").GetData(fSetSDP);
  topB.GetChild("setTimeFit").GetData(fSetTimeFit);
  topB.GetChild("setSdGeometry").GetData(fSetSdGeometry);

  // -------- do some output ------------
  // info output
  ostringstream info;
  info << " Version: "
       << GetVersionInfo(VModule::eRevisionNumber) << "\n"
    " Parameters:\n";
  info << "              set SDP: " << fSetSDP << "\n"
    "         set time fit: " << fSetTimeFit << "\n"
    "      set SD geometry: " << fSetSdGeometry << "\n";
  INFO(info);

  return eSuccess;
}



// --------------------------------------------------------------------------
VModule::ResultFlag
UseMcGeometry::Run(evt::Event& event)
{
  if (!event.HasSimShower()) {
    ERROR("No sim info. Cannot use MC geometry!");
    return eBreakLoop;
  }
  const evt::ShowerSimData& simShower = event.GetSimShower();

  if (!simShower.HasDirection()) {
    ERROR("Sim info has no direction");
    return eBreakLoop;
  }

  if (!simShower.HasPosition()) {
    ERROR("Sim info has no position");
    return eBreakLoop;
  }

  if (fSetSDP)
    SetSDP(event);

  if (fSetTimeFit)
    SetTimeFit(event);

  if (fSetSdGeometry)
    SetSdGeometry(event);

  return eSuccess;
}



// -------------------------------------------------------------------------
void
UseMcGeometry::SetSDP(evt::Event& event)
{
  if (!event.HasFEvent())
    return;

  const CoordinateSystemPtr& siteCS =
    det::Detector::GetInstance().GetSiteCoordinateSystem();

  FEvent& fevent = event.GetFEvent();
  evt::ShowerSimData& simShower = event.GetSimShower();

  const fdet::FDetector& detFD = det::Detector::GetInstance().GetFDetector();

  const CoordinateSystemPtr showerCS = simShower.GetShowerCoordinateSystem();
  const utl::Point showerCore(0,0,0,showerCS);
  utl::Vector showerAxisReal(0,0,-1,showerCS); // REAL SHOWER MOVEMENT
  if ((showerAxisReal.GetTheta(siteCS)>90.*deg))
    showerAxisReal *= -1.0; // reco-def (ALWAYS POINTING UPWARDS)

  ostringstream info;
  info << "   showerAxis theta: " << simShower.GetDirection().GetTheta(siteCS)/deg << " deg (sim) ";
  INFO(info);

  for (FEvent::EyeIterator eyeiter = fevent.EyesBegin(ComponentSelector::eHasData);
       eyeiter != fevent.EyesEnd(ComponentSelector::eHasData);
       ++eyeiter){

    fevt::Eye& eye = *eyeiter;

    if (!eye.HasHeader())
      continue;

    const fdet::Eye& eyeDet = detFD.GetEye(eye);
    const CoordinateSystemPtr eyeCS = eyeDet.GetEyeCoordinateSystem();
    utl::Point eyePos = eyeDet.GetPosition();

    utl::Vector showerAxis = showerAxisReal;
    if (showerAxis.GetTheta(eyeCS)>90.*deg)
      showerAxis *= -1; // reco-def

    if (!eye.HasRecData())
      eye.MakeRecData();

    const utl::Vector& eyeToCore = showerCore-eyePos;
    utl::Vector sdp = cross(showerAxis, eyeToCore);
    sdp.Normalize();

    EyeRecData& eyerecdata = eye.GetRecData();
    eyerecdata.SetSDPFitChiSquare(1.,1);  // reset
    eyerecdata.SetSDP(sdp);
    eyerecdata.SetSDPThetaError(0);
    eyerecdata.SetSDPPhiError(0);

    ostringstream info2;
    info2 << "          sdp theta: " << sdp.GetTheta(eyeCS)/deg
          << "                phi: " << sdp.GetPhi(eyeCS)/deg;
    INFO(info2);

    //--------------------------------------------
    // Fill Telescope geometry parameters to deal with compound eyes
    for (fevt::Eye::TelescopeIterator telIt = eye.TelescopesBegin(ComponentSelector::eHasData);
         telIt != eye.TelescopesEnd(ComponentSelector::eHasData); ++telIt)
      {
        const fdet::Telescope& detTel  = detFD.GetTelescope(*telIt);
        const utl::Point& telPos       = detTel.GetPosition();
        // 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);
        const utl::CoordinateSystemPtr telCS = detTel.GetTelescopeCoordinateSystem();

        if (!telIt->HasRecData())
          telIt->MakeRecData();
        fevt::TelescopeRecData& telRecData = telIt->GetRecData();

        const Vector telToCore = showerCore - telPos;
        //const double distCoreTel = telToCore.GetMag();

        utl::Vector showerAxis = showerAxisReal;
        if (showerAxis.GetTheta(localTelCS) > 90.*deg)
          showerAxis *= -1; // reco-def

        utl::Vector sdp = cross(showerAxis, telToCore);
        sdp.Normalize();

        telRecData.SetSDPFitChiSquare(1.,1);  // reset
        telRecData.SetSDP(sdp);
        telRecData.SetSDPThetaError(0);
        telRecData.SetSDPPhiError(0);
      }

    // loop over pixels and fill SDP-fit-pixels
    int nPix = 0;
    for (EyeRecData::PixelIterator pixeliter = eyerecdata.PulsedPixelsBegin();
         pixeliter != eyerecdata.PulsedPixelsEnd();
         ++pixeliter){

      const utl::Vector sdp = eye.GetTelescope(pixeliter->GetTelescopeId()).GetRecData().GetSDP();

      const fdet::Pixel& detpixel =
        det::Detector::GetInstance().GetFDetector().GetPixel(*pixeliter);

      const utl::Vector& dir = detpixel.GetDirection();
      const double angle = fabs(kPi/2 - fabs(Angle(dir, sdp)));

      if (angle < 2*degree) {
        eyerecdata.AddSDPPixel(*pixeliter);
        nPix++;
      } // for pixel
    }

    ostringstream info3;
    info3 << "         sdp pixels: " << nPix;
    INFO(info3);

  }// for eyeiter

}



// ------------------------------------------------------------------------------
void
UseMcGeometry::SetTimeFit(evt::Event& event)
{
  const CoordinateSystemPtr siteCS = det::Detector::GetInstance().GetSiteCoordinateSystem();

  if (!event.HasFEvent())
    return;

  FEvent&  fevent = event.GetFEvent();
  const evt::ShowerSimData& simShower = event.GetSimShower();
  const CoordinateSystemPtr showerCS = simShower.GetShowerCoordinateSystem();
  const utl::Point showerCore(0,0,0,showerCS);
  
  const TimeStamp& timeAtCore = simShower.GetTimeStamp();
  const utl::Vector showerAxisReal(0,0,-1,showerCS);

  ostringstream info;
  info << "      sim Axis (Theta,Phi) (siteCS) : "
       << showerAxisReal.GetTheta(siteCS)/deg << ", "
       << showerAxisReal.GetPhi(siteCS)/deg << " [deg]"
       << endl;

  info << "      sim Core (x,y,z)     (siteCS) : " << showerCore.GetX(siteCS)/m
       << ", " << showerCore.GetY(siteCS)/m
       << ", " << showerCore.GetZ(siteCS)/m << " [m]"
       << endl;

  const CoordinateSystemPtr localCS = simShower.GetLocalCoordinateSystem();
  const double zenith = (-simShower.GetDirection()).GetTheta(localCS);
  const double azimuth = (-simShower.GetDirection()).GetPhi(localCS);

  info << "                    simulated theta : " << zenith/deg
       << "  phi : " << azimuth/deg << " [deg]";
  INFO(info);


  for (fevt::FEvent::EyeIterator eyeiter= fevent.EyesBegin(ComponentSelector::eHasData);
       eyeiter != fevent.EyesEnd(ComponentSelector::eHasData);
       ++eyeiter) {

    fevt::Eye& eye = *eyeiter;

    if (!eye.HasHeader() || !eye.HasRecData())
      continue;

    fevt::EyeRecData& eyerecdata = eye.GetRecData();
    const utl::Vector& sdp = eyerecdata.GetSDP();

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

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

    utl::Vector showerAxis = showerAxisReal;
    bool upward = true;
    if (showerAxis.GetTheta(eyeCS)>90.*deg) {
      showerAxis *= -1; // reco-def
      upward = false;
    }

    if (upward) {
      INFO("         upward: YES");
    }
    else {
      INFO("         upward: NO");
    }

    const Vector vertical(0, 0, 1, eyeCS);
    Vector horizontalWithinSDP = cross(vertical, sdp); // horizontal
    horizontalWithinSDP.Normalize();
    if (horizontalWithinSDP.GetPhi(eyeCS) < 0 ||
        horizontalWithinSDP.GetPhi(eyeCS) > kPi) {
      //sdp *= -1;
      horizontalWithinSDP *= -1;
    }

    //VN: deal with back-eye core
    Vector eyeToCoreInEye = eyeToCore-eyeToCore.GetZ(eyeCS)/showerAxis.GetZ(eyeCS)*showerAxis;
    if (eyeToCoreInEye.GetPhi(eyeCS) < 0 ||
        eyeToCoreInEye.GetPhi(eyeCS) > kPi) {
      horizontalWithinSDP *= -1;
    }

    const double chi0 = Angle(horizontalWithinSDP, showerAxis) -
      (upward ? 180.*deg : 0.*deg);

    const double Rp = cross(showerAxis, eyeToCore).GetMag() /
       showerAxis.GetMag() * (upward ? -1 : 1);

    /*
      double delta = (showerCore.GetZ(eyeCS) - 0) / showerAxis.GetZ(eyeCS);
      utl::Point point = showerCore + showerAxis*delta;
      double distance = (point-eyePos).GetMag ();
      double Rp = sin(chi0) * distance;
    */
    /*
      utl::TimeInterval deltaT = (cos(chi0) * distance / kSpeedOfLight);
      utl::TimeStamp time = timeAtCore + deltaT;
      double T0 = (time-eyeTriggerTime).GetInterval();
    */

    const double distCoreEye = eyeToCore.GetMag();
    const double beta = Angle(-showerAxisReal, eyeToCore);
    const double T0 = ((timeAtCore + TimeInterval(distCoreEye/kSpeedOfLight * cos(beta)))
                       - eyeTriggerTime).GetInterval();

    eyerecdata.SetChiZero(chi0, 1.e-5*chi0);
    eyerecdata.SetTZero(T0, 1.e-5*T0);
    eyerecdata.SetRp(Rp, 1.e-5*Rp);

    ostringstream info;
    info << "         chi0: " << chi0/deg << " deg "
         << "           Rp: " << Rp/m << " m"
         << "           T0: " << T0  << " ns";
    INFO(info);

    //--------------------------------------------
    // Fill Telescope geometry parameters to deal with compound eyes
    for (fevt::Eye::TelescopeIterator telIt = eye.TelescopesBegin(ComponentSelector::eHasData);
         telIt != eye.TelescopesEnd(ComponentSelector::eHasData); ++telIt)
      {
        const fdet::Telescope& detTel  = detFD.GetTelescope(*telIt);
        const utl::Point& telPos       = detTel.GetPosition();
        // 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);
        const utl::CoordinateSystemPtr telCS = detTel.GetTelescopeCoordinateSystem();

        if (!telIt->HasRecData())
          telIt->MakeRecData();
        fevt::TelescopeRecData& telRecData = telIt->GetRecData();

        const Vector telToCore = showerCore - telPos;
        const double distCoreTel = telToCore.GetMag();
        const double beta = Angle(-showerAxisReal, telToCore);
        const double T0 = ((timeAtCore + TimeInterval(distCoreTel/kSpeedOfLight * cos(beta)))
                           - eyeTriggerTime).GetInterval();

        utl::Vector showerAxis = showerAxisReal;
        bool upward = true;
        if (showerAxis.GetTheta(localTelCS) > 90.*deg) {
          showerAxis *= -1; // reco-def
          upward = false;
        }

        utl::Vector sdp = cross(showerAxis, telToCore);
        sdp.Normalize();

        const Vector vertical(0, 0, 1, localTelCS);
        Vector horizontalWithinSDP = cross(vertical, sdp);
        horizontalWithinSDP.Normalize();

        // x-axis pointing towards optical axis, z-axis still vertical
        const CoordinateSystemPtr auxCS =
          CoordinateSystem::RotationZ(detTel.GetAxis().GetPhi(localTelCS), localTelCS);
        const utl::Vector horizontalProjOptAxis(1, 0, 0, auxCS);

        // Make sure that the horiz-in-sdp vector points to the front side of the
        // telescope. This is what the equivalent code above does for the eye.
        // Like the original, this is likely to blow up for core-behind-telescope cases?
        const double angle = Angle(horizontalProjOptAxis, horizontalWithinSDP);
        if (angle > kPiOnTwo)
          horizontalWithinSDP *= -1;

        //VN: deal with back-eye core
        Vector telToCoreInEye = telToCore-telToCore.GetZ(auxCS)/showerAxis.GetZ(auxCS)*showerAxis;
        const double angle2 = Angle(horizontalProjOptAxis, telToCoreInEye);
        if (angle2 > kPiOnTwo)
          horizontalWithinSDP *= -1;

        const double chi0 = Angle(horizontalWithinSDP, showerAxis)
          - (upward ? 180.*deg : 0.*deg);

        const double Rp = cross(showerAxis, telToCore).GetMag()
          / showerAxis.GetMag()
          * (upward ? -1 : 1);

        telRecData.SetChiZero(chi0, 1.e-5*chi0);
        telRecData.SetTZero(T0, 1.e-5*T0);
        telRecData.SetRp(Rp, 1.e-5*Rp);

        telRecData.SetTimeFitChiSquare(1,1); // FIXME find more elegant solution
        //telRecData.SetTimeFitCorrelations(rRpT0,rRpChi0,rChi0T0);
      }

    if (!eyerecdata.HasFRecShower())
      eyerecdata.MakeFRecShower();
    ShowerFRecData& recshower = eyerecdata.GetFRecShower();

    recshower.SetAxis(showerAxis*(upward?-1:1));

    // M.U.: FD core: intercept axis-horizontal in eye CS
    const Vector core_eye_vec = Rp / sin(kPi - chi0) * horizontalWithinSDP;
    const Point coreInEye = eyePos + core_eye_vec;
    recshower.SetCorePosition(coreInEye);

    // V.N.: Correction for timeAtCore
    const Vector coreToCoreInEye = coreInEye-showerCore;
    double timeCorrection = -coreToCoreInEye.GetMag()/kSpeedOfLight;
    if (coreToCoreInEye.GetTheta(eyeCS) > 90*deg)
      timeCorrection *= -1;
    if (upward)
      timeCorrection *= -1;
    const TimeStamp& timeAtCoreInEye = timeAtCore + timeCorrection;
    recshower.SetCoreTime(timeAtCoreInEye, 0.);

    // loop over pixels and fill SDP-fit-pixels to time-fit-pixels
    unsigned int ndof=0;
    for (EyeRecData::PixelIterator pixeliter = eyerecdata.SDPPixelsBegin ();
         pixeliter != eyerecdata.SDPPixelsEnd ();
         ++pixeliter){

      const fdet::Pixel& detpixel = detFD.GetPixel(*pixeliter);
      const fdet::Channel& detChannel = detFD.GetChannel(*pixeliter);

      // time correction for different eyes
      const int telid = pixeliter->GetTelescopeId();
      const unsigned int timeoffset = eye.GetTelescope(telid).GetTimeOffset();

      const double timebinsize = detChannel.GetFADCBinSize();
      const double t_i = pixeliter->GetRecData().GetCentroid()
        * timebinsize   // convert to nanosecond
        + timeoffset;   // add offset

      const double t_i_Err = pixeliter->GetRecData().GetCentroidError()
        * timebinsize;

      const Vector& pixeldir = detpixel.GetDirection();
      //pixeldir.TransformTo(CS);
      //sdp.TransformTo (CS);
      //fSDP.Normalize();

      //Vector vertical (0,0,1,eyeCS,Vector::kCartesian);
      //Vector horizontalWithinSDP = cross (sdp, vertical);
      //horizontalWithinSDP.Normalize();

      // check that it is the one looking outward
      //if (Angle(horizontalWithinSDP,pixeldir) >kPi/2)
      //        horizontalWithinSDP*=-1;

      const fdet::Telescope& detTel  = detFD.GetEye(eye).GetTelescope(telid);
      const utl::Point& telPos       = detTel.GetPosition();
      // 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);

      const Vector telToCore = showerCore - telPos;

      const utl::Vector sdp = eye.GetTelescope(telid).GetRecData().GetSDP();
      const Vector vertical(0, 0, 1, localTelCS);
      Vector horizontalWithinSDP = cross(vertical, sdp);
      horizontalWithinSDP.Normalize();

      const CoordinateSystemPtr auxCS =
          CoordinateSystem::RotationZ(detTel.GetAxis().GetPhi(localTelCS), localTelCS);
      const utl::Vector horizontalProjOptAxis(1, 0, 0, auxCS);

      // Make sure that the horiz-in-sdp vector points to the front side of the
      // telescope. This is what the equivalent code above does for the eye.
      // Like the original, this is likely to blow up for core-behind-telescope cases?
      const double angle = Angle(horizontalProjOptAxis, horizontalWithinSDP);
      if (angle > kPiOnTwo)
        horizontalWithinSDP *= -1;

      //VN: deal with back-eye core
      Vector telToCoreInEye = telToCore-telToCore.GetZ(auxCS)/showerAxis.GetZ(auxCS)*showerAxis;
      const double angle2 = Angle(horizontalProjOptAxis, telToCoreInEye);
      if (angle2 > kPiOnTwo)
        horizontalWithinSDP *= -1;

      // subtract component parallel to sdp
      Vector chi_i_vector = pixeldir - (sdp * pixeldir) *  sdp;
      //chi_i_vector.TransformTo(eyeCS);
      chi_i_vector.Normalize();

      const double chi_i = Angle(chi_i_vector, horizontalWithinSDP);

      pixeliter->GetRecData().SetChi_i(chi_i);
      pixeliter->GetRecData().SetT_i(t_i,t_i_Err);

      // check residual
      double chi0 = eye.GetTelescope(telid).GetRecData().GetChiZero();
      double T0 = eye.GetTelescope(telid).GetRecData().GetTZero();
      double Rp = eye.GetTelescope(telid).GetRecData().GetRp();
      const double t_expected =
        T0 + Rp/kSpeedOfLight * std::tan((chi0 - chi_i) / 2);

      ostringstream info;
      info << " chi_i=" << chi_i/deg
           << " t_i=" << t_i
           << " chi0=" << chi0
           << " t0=" << T0
           << " Rp=" << Rp
           << " t_exp=" << t_expected;
      INFO(info);


      if (fabs(t_i - t_expected) > 400)
        continue;

      if (this->IsBadPixel(*pixeliter,eye))
        continue;

      // if pixel is distant from SDP, do not take it into account
      //double distanceToSdp = abs(kPi/2 -  std::acos(fSDP * pixeldir)) ;
      //if (distanceToSdp > fTolerance) continue;

      eyerecdata.AddTimeFitPixel (*pixeliter);
      ndof++;
    }

    if (ndof)
      eyerecdata.SetTimeFitChiSquare(ndof, ndof);
    else
      eyerecdata.SetTimeFitChiSquare(1, 1);

    ostringstream info3;
    info3 << "         ndof: " << ndof;
    INFO(info3);

  } // for eyeiter
}



// ------------------------------------------------------------------
void
UseMcGeometry::SetSdGeometry(evt::Event& event)
{
  if (!event.HasSEvent())
    return;

  const CoordinateSystemPtr& siteCS =
    det::Detector::GetInstance().GetSiteCoordinateSystem();

  evt::ShowerSimData& simShower = event.GetSimShower();
  const CoordinateSystemPtr showerCS = simShower.GetShowerCoordinateSystem();
  const utl::Point showerCore(0,0,0,showerCS);
  utl::Vector showerAxisReal(0,0,-1,showerCS); // REAL SHOWER MOVEMENT
  if ((showerAxisReal.GetTheta(siteCS)>90.*deg))
    showerAxisReal *= -1.0; // reco-def (ALWAYS POINTING UPWARDS)
  const TimeStamp& timeAtCore = simShower.GetTimeStamp();


  if (!event.HasRecShower())
    event.MakeRecShower();
  if (!event.GetRecShower().HasSRecShower())
    event.GetRecShower().MakeSRecShower();
  ShowerSRecData& showersrec = event.GetRecShower().GetSRecShower();

  showersrec.SetAxis(showerAxisReal);

  //showersrec.SetAngleErrors(axis.GetCoordinates(fgLocalCS),
  //Vector::Triple(fit->sigmaU2, fit->sigmaUV, fit->sigmaV2));

  showersrec.SetCurvature(0, 0);
  showersrec.SetBarycenter(showerCore);
  showersrec.SetCorePosition(showerCore);
  showersrec.SetCoreTime(timeAtCore, TimeInterval(0));

  /*
    double timeMean3d;
    double timeSpread3d;
    double chi2;
    CalculateTimeResidual3D(axis, fit->ct0, timeMean3d, timeSpread3d, chi2);
    showersrec.SetTimeResidualMean(timeMean3d);
    showersrec.SetTimeResidualSpread(timeSpread3d);
  */

  showersrec.SetAngleChi2(0, 0);
  showersrec.SetLDFRecStage(0);
}


VModule::ResultFlag
UseMcGeometry::Finish()
{
  return eSuccess;
}



// --------------------------------------------------------------------------
bool
UseMcGeometry::IsBadPixel(const fevt::Pixel& pixel, const fevt::Eye& eye)
{
  const fdet::FDetector& detFD = det::Detector::GetInstance().GetFDetector();
  const EyeRecData& recdata = eye.GetRecData();

  // check S/N limit

  double signal = pixel.GetRecData().GetTotalCharge ();
  double rms = pixel.GetRecData().GetRMS();
  int i1 = pixel.GetRecData().GetPulseStart();
  int i2 = pixel.GetRecData().GetPulseStop();
  double nTim = (double) (i2-i1+1);
  static const double kPhoton2Pe = 0.128; // M.U. same as in FDSim
  double noise = nTim*rms*rms+signal/kPhoton2Pe;
  double sigToNoise = signal/sqrt(noise);
  static double minSigNoise = 1;
  if(sigToNoise<minSigNoise)
    return true;

  // if error is <=0. centroid finding failed
  const double t_i_Err = pixel.GetRecData().GetCentroidError();
  if (t_i_Err<=0.)
    return true;

  // check time isolation
  fevt::EyeRecData::ConstPixelIterator iter;
  for (iter = recdata.PulsedPixelsBegin();
       iter != recdata.PulsedPixelsEnd();
       ++iter){

    if (*iter==pixel)
      continue;

    const double fadcBinSize = detFD.GetChannel(*iter).GetFADCBinSize();
    const double timediff =
      std::fabs(pixel.GetRecData().GetCentroid() -
                iter->GetRecData().GetCentroid()) * fadcBinSize;

    static double isolationLimit = 10000*ns;
    if (timediff < isolationLimit )
      return false;

  } // for

  return true;
}


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