FdAxisFinder.cc

Go to the documentation of this file.

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

#include <utl/Reader.h>
#include <utl/ErrorLogger.h>
#include <utl/UTMPoint.h>

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

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

#include <det/Detector.h>

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

#include <utl/MathConstants.h>
#include <utl/PhysicalConstants.h>
#include <utl/Transformation.h>

//========================================================
// added for fitting the axis fixing the core to the clf

#include <utl/ReferenceEllipsoid.h>
#include <fwk/LocalCoordinateSystem.h>
#include <evt/ShowerRecData.h>

//========================================
// sd part for identifying Celeste or Ramiro station
// in order to identify CLF and XLF events.

#include <evt/ShowerSRecData.h>
#include <sevt/SEvent.h>
#include <sevt/Station.h>
#include <sevt/StationTriggerData.h>
#include <sevt/StationRecData.h>
#include <sevt/EventTrigger.h>
#include <sdet/SDetector.h>
#include <sdet/Station.h>

//========
//  to extract the GPS time of the event

#include <evt/Header.h>

//==================



#include <TMinuit.h>

#include "FdAxisFinder.h"

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


fevt::Eye * FdAxisFinder::fCurEye =0;
double FdAxisFinder::fChi2 =0;
unsigned int FdAxisFinder::fNDof=0;

utl::Vector FdAxisFinder::fSDP = utl::Vector();

int FdAxisFinder::fDebuging;

//------------------------------------
// For reconstructing lasers
double FdAxisFinder::fRp =0;
bool FdAxisFinder::fIsCeleste = false;
bool FdAxisFinder::fIsRamiro = false;
//------------------------------------


VModule::ResultFlag FdAxisFinder::Init(){

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

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

  topB.GetChild("ForDebug").GetData(fDebuging);
  topB.GetChild("tolerance").GetData(fTolerance);
  topB.GetChild("minPixels").GetData(fMinPixels);
  topB.GetChild("sigmasLimit").GetData(fSigmasLimit);
  topB.GetChild("SelectAndRecLasers").GetData(fSelectAndRecLasers);

  int timeFitModelSwitch;
  topB.GetChild("TimeFitModel").GetData(timeFitModelSwitch);

  fdet::FTimeFitModel& timeFitModel = fdet::FTimeFitModel::GetInstance();
  switch (timeFitModelSwitch) {
  case 0:
    timeFitModel.SetModel(fdet::FTimeFitModel::eVacuumAtm);
    break;
  case 1:
    timeFitModel.SetModel(fdet::FTimeFitModel::eRealisticAtm);
    break;
  default:
    ERROR("Invalid input parameter for TimeFitModel");
    return eFailure;
  }

  fMinimumTimeExtension = 500 * ns;
  fIsBlob = false;

  // info output
  ostringstream info;
  info << "Parameters:\n"
          "            tolerance: " << fTolerance / degree <<" deg.\n"
          "       min no. pixels: " << fMinPixels << "\n"
          "  t_i deviation limit: " << fSigmasLimit <<" sigmas";
  INFO(info);

  return eSuccess;
}


VModule::ResultFlag
FdAxisFinder::Run(evt::Event& event)
{
  if (fSelectAndRecLasers && !CheckForLasers(event))
    return eContinueLoop;

  fMinuitFailed = false;
  fIsBlob = false;
  if (!event.HasFEvent())
    return eSuccess;

  fevt::FEvent& fevent = event.GetFEvent();
  fevt::FEvent::EyeIterator eyeiter;

  unsigned int neyes = 0; // number of eyes with enough pixels
  for (eyeiter = fevent.EyesBegin(ComponentSelector::eHasData);
       eyeiter != fevent.EyesEnd(ComponentSelector::eHasData);
       ++eyeiter){

    fevt::Eye& eye = *eyeiter;

    if (!eye.HasRecData())
      continue;
    fevt::EyeRecData& eyerecdata = eye.GetRecData();

    fSDP = eyerecdata.GetSDP();

    fCurEye = &eye; // acrobacy for TMinuit

    this->FillPoints(eye);

    // this->FillPoints(eye) is needed for virtual eyes which by definition have no SDP!
    if (!fSDP.GetMag2()) {
      WARNING("This eye has no valid SDP!");
      continue;
    }

    // skip virtual eyes
    if (det::Detector::GetInstance().GetFDetector().GetEye(eye).IsVirtual())
      continue;

    if(eyerecdata.GetNTimeFitPixels() < fMinPixels)
      continue;
    if(fIsBlob== true) {
      cout <<" All pixels in this Eye were triggered within "<< fMinimumTimeExtension<<" ns\n";
      cout <<" This event may be a Blob in the mirror \n";
    }


    AxisFirstGuess(eye);
    //on the basis of time difference from calculated using first guess in function, remove outliers
    RemoveOutliers(eye);

    if (fSelectAndRecLasers)
      FindAxisLaser(eye);
    else
      FindAxis(eye);

    // Now we use the fitted Axis
    // to check whether we need to readmit
    // any pixel that was rejected when
    // using the axis first guess.

    unsigned int nopixels = eyerecdata.GetNTimeFitPixels();
    while ( ReadmitPixel(eye)) {
      if (fSelectAndRecLasers)
        FindAxisLaser(eye);
      else
        FindAxis(eye);  // Fits again the axis

      if (fDebuging){
        cout<<"\n";
        cout <<eyerecdata.GetNTimeFitPixels()-nopixels<<" pixels were re-admitted.\n";
        cout<<"\n";
      }
    }

    //rejecting pixels, one by one, if they are more than
    //fSigmasLimit away
    while ( RejectPixel(eye)) {
      if (fSelectAndRecLasers)
        FindAxisLaser(eye);
      else
        FindAxis(eye);  // Fits again the axis


    }

    //--------------------------------------
    // Defining the range for chi0 [-180,180]
    // for down goin events chi0 [0,180]
    // for up goin events chi0 [-180,0]
    if ( abs (fChi0) > kPi ) {
      fChi0 -= int ( fChi0 / (2* kPi) )* 2*kPi;
      if ( fChi0 < 0 ) fChi0 += 2* kPi;
      if ( fChi0 > kPi ) fChi0 -=  2*kPi;
    }

    //-----------------------------------------

    //--------------------------------------------
    eyerecdata.SetChiZero(fChi0,fChi0Error);
    eyerecdata.SetTZero(fT0,fT0Error);
    eyerecdata.SetRp(fRp,fRpError);
    eyerecdata.SetTimeFitChiSquare(fChi2,fNDof);

    // correlation coefficiens

    double rRpT0   = fCovariance[2][1]/fRpError/fT0Error;
    double rRpChi0 = fCovariance[2][0]/fRpError/fChi0Error;
    double rChi0T0 = fCovariance[0][1]/fChi0Error/fT0Error;

    eyerecdata.SetTimeFitCorrelations(rRpT0,rRpChi0,rChi0T0);

    //--------------------------------------------
    // Also fill the TelescopeRecData with the same values
    // since this module isn't able to deal with compound eyes
    // that consist of telescopes at multiple different positions

    for (fevt::Eye::TelescopeIterator telIt = eye.TelescopesBegin(ComponentSelector::eHasData);
         telIt != eye.TelescopesEnd(ComponentSelector::eHasData); ++telIt)
    {
      if (!telIt->HasRecData())
        telIt->MakeRecData();
      fevt::TelescopeRecData& telRecData = telIt->GetRecData();
      telRecData.SetChiZero(fChi0,fChi0Error);
      telRecData.SetTZero(fT0,fT0Error);
      telRecData.SetRp(fRp,fRpError);
      telRecData.SetTimeFitChiSquare(fChi2,fNDof);
      telRecData.SetTimeFitCorrelations(rRpT0,rRpChi0,rChi0T0);
    }

    //--------------------------------------------
    // calculate axis parameters.

    CoordinateSystemPtr eyeCS =
      det::Detector::GetInstance().GetFDetector().GetEye(eye).GetEyeCoordinateSystem();

    Vector sdp = eyerecdata.GetSDP();
    //if (sdp.GetPhi(eyeCS) < 0  ) sdp *=-1;

    Vector  vertical(0,0,1,eyeCS);
    Vector horizontalInSDP = cross(sdp, vertical); // horizontal
    horizontalInSDP.Normalize();

    Transformation rot (Transformation::Rotation(-fChi0, sdp, eyeCS));

    Vector axis (rot* horizontalInSDP);             // apply rotation
    axis.Normalize();

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

    {
      ostringstream info;
      info << "Axis(Theta,Phi)(siteCS): " << axis.GetTheta(siteCS)/degree
           << "," << axis.GetPhi(siteCS)/degree << " deg" <<endl;
      INFO(info);
    }

    Vector core_eye_vec = fRp / sin( kPi - fChi0) * horizontalInSDP;

    Point eyeposition =
      det::Detector::GetInstance().GetFDetector().GetEye(eye).GetPosition();

    Point core = eyeposition + core_eye_vec;

    if (fDebuging == 1){
      cout << "Core (x,y)(eyeCS) " << core.GetX(eyeCS)<<","
           <<  core.GetY(eyeCS)<< " [m]" << endl;

      cout << "Core (x,y)(siteCS) " << core.GetX(siteCS)<<","
           <<  core.GetY(siteCS)<< " [m]" << endl;
    }

    if (!eyerecdata.HasFRecShower())
      eyerecdata.MakeFRecShower();

    ShowerFRecData& recshower = eyerecdata.GetFRecShower();
    recshower.SetAxis(axis);
    recshower.SetCorePosition(core);
    recshower.SetCoreTime(fCurEye->GetHeader().GetTimeStamp() -
                          TimeInterval(1000*100*ns), fRp, fChi0, fT0);

    ++neyes;

    {
      ostringstream info;
      info << "T_0=" << fT0/ns  << " +/- " << fT0Error/ns << " ns,  "
           << "Chi_O="  << fChi0/degree << " +/- " << fChi0Error/degree << " deg, "
           << "Rp   "  << fRp/meter   << " +/- " << fRpError/meter << " m";
      INFO(info);
    }

  }// for eyeiter

  if (!neyes)
    return eContinueLoop; //no eyes worth reconstructing, next event

  return eSuccess;
}


void
FdAxisFinder::FindAxis(fevt::Eye& /*eye*/)
{
  const int npar = 3; // number of parameters

  TMinuit theMinuit(npar); // 3 parameter fit

  theMinuit.SetPrintLevel(-1); // minuit verbosity
  theMinuit.SetFCN(FdAxisFinder::MinuitFitFunc);

  double arglist[npar];  // arguments to pass to the Minuit interpreter
  int ierflag = 0;

  arglist[0] = 1; // [up]
  arglist[1] = 1; //  unused
  theMinuit.mnexcm("SET ERR", arglist, 1, ierflag);

  if (ierflag)
    ERROR("Minuit SET ERR failed");

  static double vstart[npar]; // inital point
  static double step[npar];   // step size

  // something better has to be found for starting points
  vstart[0] = fChi0First;
  vstart[1] = fT0First;
  vstart[2] = fRpFirst;

  step[0] = 0.001;
  step[1] = 10;
  step[2] = 10;

  theMinuit.mnparm(0, "Chi0", vstart[0], step[0], 0, 0, ierflag);
  theMinuit.mnparm(1, "T0", vstart[1], step[1], 0, 0, ierflag);
  theMinuit.mnparm(2, "Rp", vstart[2], step[2], 0, 0, ierflag);

  arglist[0] = 3000;  // [maxcalls]
  arglist[1] = 1;    // [tolerance]
  theMinuit.mnexcm("MIGRAD", arglist, 2, ierflag);

  if (ierflag) {
    ERROR("Minuit MIGRAD failed");
    fMinuitFailed = true;
  }

  theMinuit.mnexcm("HESSE", arglist, 2, ierflag);

  // Check the status of the covariance matrix
  // if bad use MINOS
  double fmin, fedm, errdef; // unused dummies
  int npari, nparx; //unused dummies
  int istat = 0;
  theMinuit.mnstat(fmin, fedm, errdef, npari, nparx, istat);
  if (istat < 3) { // no accurate covariance matrix
    theMinuit.mnexcm("MINOS", arglist, 2, ierflag);
    theMinuit.mnstat(fmin, fedm, errdef, npari, nparx, istat);
    if (istat < 3)
      WARNING("Neither HESSE nor MINOS found accurate covariance matrix");
  }
  theMinuit.GetParameter(0, fChi0, fChi0Error);
  theMinuit.GetParameter(1, fT0, fT0Error);
  theMinuit.GetParameter(2, fRp, fRpError);
  theMinuit.mnemat(&fCovariance[0][0], 3);

  if (fDebuging == 1) {
    cout << "TO   "  << fT0  << " +/- " << fT0Error << "\n"
            "ChiO "  << fChi0/degree << " +/- " << fChi0Error/degree << "\n"
            "Rp   "  << fRp   << " +/- " << fRpError << endl;
  }

  fChi2 = theMinuit.fAmin;
}


void
FdAxisFinder::FillPoints(fevt::Eye& eye)
{
  double min_t_i = 10e6;
  double max_t_i = 0;

  fevt::EyeRecData& recdata = eye.GetRecData();

  const fdet::Eye& deteye =
    det::Detector::GetInstance().GetFDetector().GetEye(eye);

  // do all calculation in the eye reference system
  CoordinateSystemPtr CS = deteye.GetEyeCoordinateSystem();

  //------------------------------
  //all pixels used for SDP are converted into timeFit pixels
  //if they are within fTolerance.
  for (fevt::EyeRecData::PixelIterator pixeliter = recdata.SDPPixelsBegin();
       pixeliter != recdata.SDPPixelsEnd(); ++pixeliter) {

    fevt::Pixel& evtpixel = *pixeliter;
    const fdet::Pixel& detpixel =
      det::Detector::GetInstance().GetFDetector().GetPixel(evtpixel);
    const double AngleWithSDP = fabs(kPi/2 - Angle(fSDP, detpixel.GetDirection()));

    if (AngleWithSDP < fTolerance)
      recdata.AddTimeFitPixel(*pixeliter);

  }
  //-------------------------------

  //--------------------------------------------
  //This loop will fill the information
  //for timing fit (t_i, t_i_error and chi_i).
  //The loop is over all pixels with pulse
  for (fevt::EyeRecData::PixelIterator pixeliter = recdata.PulsedPixelsBegin();
       pixeliter != recdata.PulsedPixelsEnd(); ++pixeliter) {

    fevt::Pixel& evtpixel = *pixeliter;
    const fdet::Pixel& detpixel =
      det::Detector::GetInstance().GetFDetector().GetPixel(evtpixel);
    const fdet::Channel& detChannel =
      det::Detector::GetInstance().GetFDetector().GetChannel(evtpixel);

    double timebinsize = detChannel.GetFADCBinSize();

    // time correction for different eyes
    int telid = evtpixel.GetTelescopeId();
    unsigned int timeoffset = eye.GetTelescope(telid).GetTimeOffset();
    double t_i = evtpixel.GetRecData().GetCentroid()  ;
    t_i *= timebinsize;  // convert to nanosecond
    t_i += timeoffset;   // add offset
    const double t_i_Err = timebinsize * evtpixel.GetRecData().GetCentroidError();

    if (fDebuging == 1)
      cout << "PixId: " << evtpixel.GetId() << " t_i " << t_i << " ti err " << t_i_Err << '\n';

    // SDP is 0 for virtual eye, fill times nevertheless, but no chi
    if (fSDP.GetMag2()) {
      //Estimating horizontal Eye-Core direction
      Vector pixeldir = detpixel.GetDirection();
      pixeldir.TransformTo(CS);
      Vector vertical(0,0,1, CS, Vector::kCartesian);
      Vector horizontalWithinSDP = Cross(fSDP, vertical);
      horizontalWithinSDP /= horizontalWithinSDP.GetMag();
      // subtract component parallel to sdp
      Vector chi_i_vector = pixeldir - (fSDP * pixeldir) * fSDP;
      chi_i_vector.TransformTo(CS);
      //  Estimating chi_i
      chi_i_vector /= chi_i_vector.GetMag();
      const double chi_i = Angle(chi_i_vector, horizontalWithinSDP);

      evtpixel.GetRecData().SetChi_i(chi_i);

    } else
      evtpixel.GetRecData().SetChi_i(0);

    evtpixel.GetRecData().SetT_i(t_i, t_i_Err);

    if (t_i < min_t_i)
      min_t_i = t_i;
    if (t_i > max_t_i)
      max_t_i = t_i;

  } // for pixels with pulse

  if (abs(max_t_i -  min_t_i) < fMinimumTimeExtension)
    fIsBlob = true;
}


void
FdAxisFinder::MinuitFitFunc(int& /*npar*/, double* /*gin*/, double& f,
                            double *par, int /*iflag*/)
{
  fevt::EyeRecData& recdata = fCurEye->GetRecData();

  double chisq = 0;
  unsigned int ndof = 0;

  const double& chi_0 = par[0];
  const double& t_0 = par[1];
  const double& r_p = par[2];

  const fdet::FTimeFitModel& timeFitModel = fdet::FTimeFitModel::GetInstance();
  CoordinateSystemPtr eyeCS =
    det::Detector::GetInstance().GetFDetector().GetEye(*fCurEye).GetEyeCoordinateSystem();
  const double sdpTheta = fSDP.GetTheta(eyeCS);
  const int eyeID = fCurEye->GetId();

  for (EyeRecData::PixelIterator iter = recdata.TimeFitPixelsBegin();
       iter != recdata.TimeFitPixelsEnd(); ++iter) {

    const PixelRecData& pixrec = iter->GetRecData();
    const int telID = iter->GetTelescopeId();

    const double t_i = pixrec.GetT_i().GetNanoSecond();
    const double t_i_Err = pixrec.GetT_iError().GetInterval();
    const double chi_i = pixrec.GetChi_i();

    const double t_expected = timeFitModel.GetTimeAtAperture(t_0, r_p, chi_0, chi_i, sdpTheta, eyeID, telID);
    const double tmp = t_i - t_expected;

    chisq += tmp*tmp / (t_i_Err * t_i_Err);
    ++ndof;
  }

  f = chisq;

  fNDof = ndof -3;
}


void
FdAxisFinder::MinuitFitFuncLaser(int& /*npar*/, double* /*gin*/, double& f,
                                 double *par, int /*iflag*/)
{
  fevt::EyeRecData& recdata = fCurEye->GetRecData();

  double chisq = 0;
  unsigned int ndof = 0;

  const double& chi_0 = par[0];
  const double& t_0 = par[1];
  double& r_p = fRp; //initializing Rp

  //===========================================================
  //Estimating Rp from the SDP and chi_0 and using the CLF
  // location as core.

  CoordinateSystemPtr CS =
    det::Detector::GetInstance().GetSiteCoordinateSystem();

  const ReferenceEllipsoid& e = ReferenceEllipsoid::GetWGS84();
  //UTMPoint UTMclf(6095769.0,469378.0,1412.0,19,'H',e);

  //--------------------------------
  Point eyeposition =
    det::Detector::GetInstance().GetFDetector().GetEye(*fCurEye).GetPosition();
  UTMPoint UTMeye(eyeposition, e);
  Point ThisEye = UTMeye.GetPoint(CS);
  //------------------------

  const CoordinateSystemPtr xlfCs = fwk::CoordinateSystemRegistry::Get("XLF");
  const CoordinateSystemPtr clfCs = fwk::CoordinateSystemRegistry::Get("CLF");

  Point clf(0,0,0, clfCs);
  const UTMPoint UTMclf(clf, e);
  clf = UTMclf.GetPoint(CS);

  Point xlf(0,0,0, xlfCs);
  const UTMPoint UTMxlf(xlf, e);
  xlf = UTMxlf.GetPoint(CS);

  CoordinateSystemPtr EyeCS = fwk::LocalCoordinateSystem::Create(ThisEye);

  Point LaserFacility;
  if (fIsCeleste)
    LaserFacility = clf;
  if (fIsRamiro)
    LaserFacility = xlf;

  Vector vertical_eye(0,0,1, EyeCS);
  Vector eye_LaserFacility = LaserFacility - ThisEye;

  //----------------------------

  Vector horizontalInSDP = Cross(fSDP, vertical_eye);
  horizontalInSDP /= horizontalInSDP.GetR(EyeCS);

  Transformation rot1(Transformation::Rotation(-chi_0, fSDP, EyeCS));
  Vector LaserOrientation1(rot1 * horizontalInSDP);
  LaserOrientation1 /= LaserOrientation1.GetR(EyeCS);

  //---------------------------
  Vector newSDP = Cross(eye_LaserFacility, LaserOrientation1);

  if (abs(newSDP.GetTheta(EyeCS) - fSDP.GetTheta(EyeCS)) < kPi/2)
    fSDP = newSDP;
  else
    fSDP = -newSDP;

  horizontalInSDP = Cross(fSDP, vertical_eye);
  horizontalInSDP /= horizontalInSDP.GetR(EyeCS);

  Transformation rot2(Transformation::Rotation(-chi_0, fSDP, EyeCS));
  Vector LaserOrientation(rot2 * horizontalInSDP);
  LaserOrientation /= LaserOrientation.GetR(EyeCS);

  const double beta = Angle(eye_LaserFacility, LaserOrientation);

  r_p = eye_LaserFacility.GetR(EyeCS) * sin(beta);

  //=====================================

  const fdet::FTimeFitModel& timeFitModel = fdet::FTimeFitModel::GetInstance();
  CoordinateSystemPtr eyeCS =
    det::Detector::GetInstance().GetFDetector().GetEye(*fCurEye).GetEyeCoordinateSystem();
  double sdpTheta = fSDP.GetTheta(eyeCS);
  int eyeID = fCurEye->GetId();

  for (EyeRecData::PixelIterator iter = recdata.TimeFitPixelsBegin();
       iter != recdata.TimeFitPixelsEnd(); ++iter) {

    const PixelRecData& pixrec = iter->GetRecData();
    const int telID = iter->GetTelescopeId();

    const double t_i = pixrec.GetT_i().GetNanoSecond();
    const double t_i_Err = pixrec.GetT_iError().GetInterval();
    const double chi_i = pixrec.GetChi_i();

    const double t_expected = timeFitModel.GetTimeAtAperture(t_0, r_p, chi_0, chi_i, sdpTheta, eyeID, telID);
    const double tmp = t_i - t_expected;

    chisq += tmp*tmp / (t_i_Err * t_i_Err);
    ++ndof;
  }// for iter

  f = chisq;

  fNDof = ndof - 2;
}


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


void
FdAxisFinder::RemoveOutliers(fevt::Eye& eye)
{
  if (!eye.HasRecData())
    return;
  CoordinateSystemPtr eyeCS =
    det::Detector::GetInstance().GetFDetector().GetEye(eye).GetEyeCoordinateSystem();
  const fdet::FTimeFitModel& timeFitModel = fdet::FTimeFitModel::GetInstance();
  int eyeID = eye.GetId();
  double sdpTheta = fSDP.GetTheta(eyeCS);

  EyeRecData& recdata = eye.GetRecData();
  EyeRecData::PixelIterator iter;

  for (iter = recdata.TimeFitPixelsBegin();
       iter != recdata.TimeFitPixelsEnd(); ) {

    const fevt::Pixel& pixrec = *iter;
    const int telID = pixrec.GetTelescopeId();

    const double t_i = pixrec.GetRecData().GetT_i().GetNanoSecond();
    const double chi_i = pixrec.GetRecData().GetChi_i();
    const double t_i_Err = pixrec.GetRecData().GetT_iError().GetInterval();

    const double t_expected = timeFitModel.GetTimeAtAperture(fT0First, fRpFirst,  fChi0First, chi_i, sdpTheta, eyeID, telID);

    const double tmp = t_i - t_expected;
    const double time_residue = tmp / (t_i_Err) ;
    if (fabs(time_residue) > 10 && recdata.GetNTimeFitPixels() > 3)
      iter = recdata.RemoveTimeFitPixel(iter);
    else
      ++iter;

  }
}


void
FdAxisFinder::AxisFirstGuess(fevt::Eye& eye)
{
  double m[3][4] = { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } };
  int npoint = 0;

  const fevt::EyeRecData& recdata = eye.GetRecData();

  //calculate chimed
  double chimin = kPi;
  double chimax = 0.;
  for (fevt::EyeRecData::ConstPixelIterator iter = recdata.TimeFitPixelsBegin();
       iter != recdata.TimeFitPixelsEnd(); ++iter) {

    const Pixel& pixel = *iter;
    const PixelRecData& pixrec = pixel.GetRecData();
    if (pixrec.GetChi_i() < chimin)
      chimin = pixrec.GetChi_i();
    if (pixrec.GetChi_i() > chimax)
      chimax = pixrec.GetChi_i();

    ++npoint;
  }
  const double chimed = 0.5 * (chimin + chimax);

  for (fevt::EyeRecData::ConstPixelIterator iter = recdata.TimeFitPixelsBegin();
       iter != recdata.TimeFitPixelsEnd(); ++iter) {

    const Pixel& pixel = *iter;
    const PixelRecData& pixrec = pixel.GetRecData();

    const double t_i = pixrec.GetT_i().GetNanoSecond();
    const double t_i_Err = pixrec.GetT_iError().GetInterval();
    const double x = pixrec.GetChi_i() - chimed;
    //const double chi_i_Err = 0;
    const double weight = sqrt(pixrec.GetTotalCharge())/(t_i_Err*t_i_Err);

    //initial matrix
    m[0][0] += weight;
    m[0][1] += weight*x;
    m[0][2] += weight*x*x*0.5;
    m[0][3] += weight*t_i;
    m[1][0] = m[0][1];
    m[1][1] += weight*x*x;
    m[1][2] += 0.5*weight*x*x*x;
    m[1][3] += weight*t_i*x;
    m[2][0] = m[1][1];
    m[2][1] += weight*x*x*x;
    m[2][2] += 0.5*weight*x*x*x*x;
    m[2][3] += weight*t_i*x*x;
  }

  //Gauss method
  for (int k = 0; k < 2; ++k)
    for (int i = k+1; i < 3; ++i) {
      const double temp = m[i][k]/m[k][k];
      for (int j = 0; j < 4; ++j)
        m[i][j] -= temp * m[k][j];
    }

  //evaluate parameters
  const double a2 = m[2][3]/m[2][2];
  const double a1 = (m[1][3] - m[1][2]*a2)/m[1][1];
  const double a0 = (m[0][3] - m[0][1]*a1 - m[0][2]*a2)/m[0][0];
  fChi0First = chimed + 2*atan(-a2/a1);
  fRpFirst = -2*kSpeedOfLight*a1*a1*a1/(a1*a1+a2*a2);
  fT0First = a0 - 2*a2*a1*a1/(a1*a1+a2*a2);

  if (fDebuging == 1)
    cout << "First Guess: Chi0 "<< fChi0First/degree << " Rp "
         << fRpFirst << " T0 " << fT0First << endl;
}


bool
FdAxisFinder::ReadmitPixel(fevt::Eye& eye)
{
  int ReadmittedPixels = 0;
  CoordinateSystemPtr eyeCS =
    det::Detector::GetInstance().GetFDetector().GetEye(eye).GetEyeCoordinateSystem();
  const fdet::FTimeFitModel& timeFitModel = fdet::FTimeFitModel::GetInstance();
  int eyeID = eye.GetId();
  double sdpTheta = fSDP.GetTheta(eyeCS);

  EyeRecData& recdata = eye.GetRecData();

  // iterate over pixels with a pulse
  for (EyeRecData::PixelIterator pixeliter = recdata.PulsedPixelsBegin();
       pixeliter != recdata.PulsedPixelsEnd(); ++pixeliter) {

    fevt::Pixel& evtpixel = *pixeliter;

    const unsigned int pixel_id = evtpixel.GetId();
    //-----------------------------
    // Checking whether this pixel was already
    // selected for Axis Time fitting.
    bool timefit_pixel = false;
    for (EyeRecData::ConstPixelIterator iter = recdata.TimeFitPixelsBegin();
         iter != recdata.TimeFitPixelsEnd(); ++iter)
      if (pixel_id == iter->GetId()) {
        timefit_pixel = true;
        break;
      }
    //-------------------------------------
    if (timefit_pixel)
      continue;

    PixelRecData& pixrec = evtpixel.GetRecData();
    const double t_i = pixrec.GetT_i().GetNanoSecond();
    const double chi_i = pixrec.GetChi_i();
    const double t_i_Err = pixrec.GetT_iError().GetInterval();

    const int telID = evtpixel.GetTelescopeId();

    double t_expected = timeFitModel.GetTimeAtAperture(fT0, fRp,  fChi0, chi_i, sdpTheta, eyeID, telID);
    double tmp = t_i - t_expected;
    double time_residue = tmp / (t_i_Err);

    const fdet::Pixel& detpixel =
      det::Detector::GetInstance().GetFDetector().GetPixel(evtpixel);
    const double AngleWithSDP = fabs(kPi/2 - Angle(fSDP, detpixel.GetDirection()));

    //-------------------------------------
    // If the the pixel wasn't selected for timing fitting,
    // and has small time residue and is close to the SDP. The pixel
    // is added for timing fit.
    if ( fabs(time_residue) < 13 && AngleWithSDP < fTolerance) {
      recdata.AddTimeFitPixel(evtpixel);
      ReadmittedPixels++;
      if (fDebuging == 1)
        cout << "readmitting pixel  " << evtpixel.GetId() << " for the Timing Fit\n";
    }

  }
  return ReadmittedPixels > 0;
}


bool
FdAxisFinder::RejectPixel(fevt::Eye& eye)
{
  int RejectedPixels = 0;
  double residue_max = 0;
  int residue_max_pixid = -1;
  int residue_max_telid = -1;

  CoordinateSystemPtr eyeCS =
    det::Detector::GetInstance().GetFDetector().GetEye(eye).GetEyeCoordinateSystem();
  const fdet::FTimeFitModel& timeFitModel = fdet::FTimeFitModel::GetInstance();
  int eyeID = eye.GetId();
  double sdpTheta = fSDP.GetTheta(eyeCS);

  EyeRecData& eyerecdata = eye.GetRecData();

  //-------------------------------------------
  // Loop to Get the pixel that has the
  // biggest timing chisq
  for (EyeRecData::PixelIterator iter = eyerecdata.TimeFitPixelsBegin();
       iter != eyerecdata.TimeFitPixelsEnd(); ++iter) {

    const fevt::Pixel& evtpixel = *iter;
    // Estimating the contribution of each pixel to the
    // time chisq (chisq_i).
    // loading the pixel information
    const int telid = evtpixel.GetTelescopeId();
    const double t_i = evtpixel.GetRecData().GetT_i().GetNanoSecond();
    const double t_i_Err = evtpixel.GetRecData().GetT_iError().GetInterval();
    const double chi_i = evtpixel.GetRecData().GetChi_i();

    const double t_expected =
      timeFitModel.GetTimeAtAperture(fT0,fRp,fChi0,chi_i,sdpTheta,eyeID,telid);

    const double temp = t_i - t_expected;
    const double time_residue = temp / t_i_Err;
    if (fDebuging == 1)
      cout << "time chisq_i: " << (temp * temp ) / (t_i_Err * t_i_Err) << " pixelid " << evtpixel.GetId()
           << " delta t: " << temp << " ti error: " << t_i_Err << " N Sigmas away: " << time_residue  << "\n";

    if ( fabs(time_residue) > residue_max) {
      residue_max = fabs(time_residue);
      residue_max_pixid = evtpixel.GetId();
      residue_max_telid = telid;
    }
  }
  //----End find pixel with biggest time chisq

  //---------------------------------------------
  //Rejecting the pixel that may have bad pulse centroid
  //----------------------------------------------
  // if the pixel with the biggest contribution to the
  // time chisq is more than fSigmasLimit away it will be removed.
  // If there are only 3 pixels or less, no pixel is removed.
  //-----------------------------------------------------


  if (residue_max > fSigmasLimit  && eyerecdata.GetNTimeFitPixels() > 3) {

    for (EyeRecData::PixelIterator iter = eyerecdata.TimeFitPixelsBegin();
         iter != eyerecdata.TimeFitPixelsEnd(); ) {

      const fevt::Pixel& evtpixel = *iter;
      // loading the pixel information
      const int telid = evtpixel.GetTelescopeId();
      const int pixid = evtpixel.GetId();
      if (telid == residue_max_telid && pixid == residue_max_pixid) {
        iter = eyerecdata.RemoveTimeFitPixel(iter);
        RejectedPixels++;
        break;
      } else
        ++iter;
    }
  }

  return RejectedPixels > 0;
}


void
FdAxisFinder::FindAxisLaser(fevt::Eye& eye)
{
  const int npar = 2; // number of parameters

  TMinuit theMinuit(npar); // 2 parameter fit

  theMinuit.SetPrintLevel(-1); // minuit verbosity
  theMinuit.SetFCN(FdAxisFinder::MinuitFitFuncLaser);

  double arglist[npar];  // arguments to pass to the Minuit interpreter
  int ierflag = 0;

  arglist[0] = 1; // [up]
  arglist[1] = 1; //  unused
  theMinuit.mnexcm("SET ERR", arglist, 1, ierflag);

  if (ierflag)
    ERROR("Minuit SET ERR failed");

  static double vstart[npar]; // inital point
  static double step[npar];   // step size

  // something better has to be found for starting points
  vstart[0] = fChi0First;
  vstart[1] = fT0First;

  step[0] = 0.001;
  step[1] = 10;

  theMinuit.mnparm(0, "Chi0", vstart[0], step[0], 0, 0, ierflag);
  theMinuit.mnparm(1, "T0", vstart[1], step[1], 0, 0, ierflag);

  arglist[0] = 3000;  // [maxcalls]
  arglist[1] = 1;    // [tolerance]
  theMinuit.mnexcm("MIGRAD", arglist, 2, ierflag);

  if (ierflag) {
    ERROR("Minuit MIGRAD failed");
    fMinuitFailed = true;
  }
  theMinuit.GetParameter(0, fChi0, fChi0Error);
  theMinuit.GetParameter(1, fT0, fT0Error);

  //=========================================
  //Estimating Rp using the SDP and fChi0

  //===========================================================
  //Estimating Rp from the SDP and chi_0 and using the CLF
  // location as core.

  CoordinateSystemPtr CS =
    det::Detector::GetInstance().GetSiteCoordinateSystem();

  const ReferenceEllipsoid &e = ReferenceEllipsoid::GetWGS84();
  //UTMPoint UTMclf(6095769.0,469378.0,1412.0,19,'H',e);

  //--------------------------------
  Point eyeposition =
    det::Detector::GetInstance().GetFDetector().GetEye(eye).GetPosition();
  UTMPoint UTMeye(eyeposition, e);
  Point ThisEye = UTMeye.GetPoint(CS);
  //------------------------

  const CoordinateSystemPtr xlfCs = fwk::CoordinateSystemRegistry::Get("XLF");
  const CoordinateSystemPtr clfCs = fwk::CoordinateSystemRegistry::Get("CLF");

  Point clf(0,0,0, clfCs);
  const UTMPoint UTMclf(clf, e);
  clf = UTMclf.GetPoint(CS);

  Point xlf(0,0,0, xlfCs);
  const UTMPoint UTMxlf(xlf, e);
  xlf = UTMxlf.GetPoint(CS);

  CoordinateSystemPtr EyeCS = fwk::LocalCoordinateSystem::Create(ThisEye);

  Point LaserFacility;
  if (fIsCeleste)
    LaserFacility = clf;
  if (fIsRamiro)
    LaserFacility = xlf;

  Vector vertical_eye(0,0,1, EyeCS);
  Vector eye_LaserFacility = LaserFacility - ThisEye;

  //----------------------------

  Vector horizontalInSDP = cross(fSDP,vertical_eye);
  horizontalInSDP /= horizontalInSDP.GetR(EyeCS);
  Transformation rot1 (Transformation::Rotation(-fChi0,fSDP, EyeCS));
  Vector LaserOrientation1 (rot1* horizontalInSDP);
  LaserOrientation1 /= LaserOrientation1.GetR(EyeCS);

  //---------------------------
  Vector newSDP = cross(eye_LaserFacility,LaserOrientation1);

  if (abs(newSDP.GetTheta(EyeCS) - fSDP.GetTheta(EyeCS)) < kPi/2)
    fSDP = newSDP;
  else
    fSDP = -newSDP;

  horizontalInSDP = cross(fSDP,vertical_eye);
  horizontalInSDP /= horizontalInSDP.GetR(EyeCS);

  Transformation rot2 (Transformation::Rotation(-fChi0,fSDP, EyeCS));
  Vector LaserOrientation (rot2* horizontalInSDP);
  LaserOrientation /= LaserOrientation.GetR(EyeCS);

  double beta = Angle(eye_LaserFacility, LaserOrientation);

  fRp = eye_LaserFacility.GetR(EyeCS) * sin(beta);

  //=====================================

  if (fDebuging == 1) {
    cout << "TO   " << fT0 << " +/- " << fT0Error << "\n"
            "ChiO " << fChi0/degree << " +/- " << fChi0Error/degree << "\n"
            "Rp   " << fRp << " +/- " << fRpError << endl;
  }

  fChi2 = theMinuit.fAmin;
}


bool
FdAxisFinder::CheckForLasers(evt::Event& event)
{
  bool ItIsLaser = false;

  //=====================================
  // Check the GPSNanoSecond
  // (0.5 Seconds is CLF)
  // (0.7 Seconds is XLF)
  const evt::Header& header = event.GetHeader();
  const TimeStamp time = header.GetTime();
  long int GPSNano = time.GetGPSNanoSecond();
  int MaxDiff = 500000;
  bool IsCLFLaser = false;
  bool IsXLFLaser = false;

  //  hybrid lasers for SD/FD time synchronization
  if (fSelectAndRecLasers == 1) {
    if (fabs(GPSNano - 500000000) < MaxDiff)
      IsCLFLaser = true;
    if (fabs(GPSNano - 700000000) < MaxDiff)
      IsXLFLaser = true;

    if (IsCLFLaser || IsXLFLaser)
      ItIsLaser = true;
    else
      return ItIsLaser;
  }

  // vertical lasers for aerosol measurements
  if (fSelectAndRecLasers == 2) {
    if (fabs(GPSNano - 250000000) < MaxDiff) {
      IsCLFLaser = true;
      fIsCeleste = true;
      ItIsLaser = true;
      return ItIsLaser;
    } else if (fabs(GPSNano - 350000000) < MaxDiff) {
      IsXLFLaser = true;
      fIsRamiro = true;
      ItIsLaser = true;
      return ItIsLaser;
    } else {
      ItIsLaser = false;
      return ItIsLaser;
    }
  }
  //=======================================================
  //======================================================================
  //  Check whether the CLF or the XLF have trigered Celeste or Ramiro stations

  if (!event.HasSEvent()) {
    cout << " don't have SEvent" << endl;
    return false;
  }
  sevt::SEvent& sevent = event.GetSEvent();
  sevt::SEvent::StationIterator sditer;

  fIsCeleste = false;
  fIsRamiro = false;

  for (sditer = sevent.StationsBegin(); sditer != sevent.StationsEnd(); ++sditer) {
    sevt::Station& station = *sditer;

    if (!sditer->HasRecData())
      sditer->MakeRecData();

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

    if (station.IsCandidate() && station.GetId() == 203)
      fIsCeleste = true;
    if (station.IsCandidate() && station.GetId() == 805)
      fIsRamiro = true;
  }

  cout << "\n"
          "==============================================\n"
          "\n"
          "\n";
  if (fIsCeleste && IsCLFLaser) {
    cout << "        This is a CLF event\n";
    ItIsLaser = true;
  }
  if (fIsRamiro && IsXLFLaser) {
    cout << "         This is a XLF event\n";
    ItIsLaser = true;
  }
  if (!(fIsCeleste && IsCLFLaser) && !(fIsRamiro && IsXLFLaser)) {
    cout << "    This is NOT a hybrid laser event\n"
            "            (skipping event)\n";
    ItIsLaser = false;
  }
  cout << "\n"
          "\n"
          "==============================================\n"
          "\n"
          "\n";

  //-------------------
  return ItIsLaser;
}