HybridGeometryFinderOG/HybridGeometryFinder.cc

Go to the documentation of this file.

#include <iostream>
#include <fstream>
#include <iomanip>
#include <algorithm>

#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/TimeStamp.h>
#include <utl/Point.h>
#include <utl/UTMPoint.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/Math.h>

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

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

#include <det/Detector.h>

#include <fdet/FDetector.h>
#include <fdet/Eye.h>
#include <fdet/Telescope.h>
#include <fdet/Channel.h>
#include <fdet/Pixel.h>
#include <fdet/FTimeFitModel.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>

#include "HybridGeometryFinder.h"

//debug
#include <TCanvas.h>
#include <TF1.h>
#include <TGraph.h>
#include <TFile.h>
#include <TMinuit.h>

using namespace std;
using namespace fwk;
using namespace utl;
using namespace fevt;
using namespace evt;
using namespace HybridGeometryFinderOG;
using namespace sevt;


fevt::Eye* HybridGeometryFinder::fCurEye = 0;
evt::Event* HybridGeometryFinder::fCurEvent = 0;

double HybridGeometryFinder::fDistanceStationSDP = 2000.;
double HybridGeometryFinder::fVerticalCurvature = 0.;
int HybridGeometryFinder::fHotStationId = 1;
int HybridGeometryFinder::fround = 1;
int HybridGeometryFinder::fTriggeredStations = 1;
double HybridGeometryFinder::fSdHybridOffset = 0;
double HybridGeometryFinder::fStationAxisDistance = 9999.9;
double HybridGeometryFinder::fsignal_max = 0;
double HybridGeometryFinder::fsignal_max_pre = 0;


int HybridGeometryFinder::fHybridCore = 0;
int HybridGeometryFinder::fDebuging = 0;
double HybridGeometryFinder::fChi2TimeFit =0;
int HybridGeometryFinder::fNDofTimeFit =0;
double HybridGeometryFinder::fNsCorrection=0;
int HybridGeometryFinder::fCeleste =0;
utl::Vector HybridGeometryFinder::fSDP = utl::Vector();

//FS (June 2007): check for Minuit convergence
bool HybridGeometryFinder::fFoundAxis=false;


VModule::ResultFlag HybridGeometryFinder::Init(){

  CentralConfig* cc = CentralConfig::GetInstance();

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


  topB.GetChild("ForDebug").GetData(fDebuging);
  topB.GetChild("DistanceStationSDP").GetData(fDistanceStationSDP);
  topB.GetChild("VerticalCurvature").GetData(fVerticalCurvature);
  topB.GetChild("tolerance").GetData(fTolerance);
  topB.GetChild("sigmasLimit").GetData(fSigmasLimit);
  topB.GetChild("UseNthStation").GetData(fUseNthStation);
  if ( fUseNthStation == 0 ) fOptimize = true;
  else  fOptimize = false;

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

  // info output
  ostringstream info;
  info << " Version: "
       << GetVersionInfo(VModule::eRevisionNumber) << "\n"
       << " Parameters:\n"
       << "\ttank distance to sdp: " << fDistanceStationSDP/m << " m\n"
       << "\t  vertical curvature: " << fVerticalCurvature*m << " m^(-1) "
       << "\ttime fit uses " << timeFitModel.GetModelName() << "\n";
  INFO(info);

  return eSuccess;
}


VModule::ResultFlag HybridGeometryFinder::Run(evt::Event& event){

  INFO("HybridGeometryFinder::Run()");

  if (!event.HasSEvent())
    return eSuccess;

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

  //sevt::SEvent& sevent = event.GetSEvent();

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

  fevt::FEvent::EyeIterator eyeiter;


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


    fevt::Eye& eye = *eyeiter;

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

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

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

    if (eyerecdata.GetNPulsedPixels() < 3) {
       cout << "This eye "<<eye.GetId()<<" has only "<<
        eyerecdata.GetNPulsedPixels()<<" pixels with pulse \n";
      continue;
    }

    if (!eyerecdata.GetSDP().GetMag2()) {
      WARNING("This eye has no SDP!");
      continue;
    }

    bool FirstGoodSt = false;
    int nHybridStations = 0;
    int UseStation = 0;
    fSDP =  eyerecdata.GetSDP();
    fChi0 = eyerecdata.GetChiZero();
    fT0 =   eyerecdata.GetTZero();
    fRp =   eyerecdata.GetRp();

  //-------------------------------
  // If reconstructing lasers a correction of +180 ns
  // is required to apply to the FD time.
  // The time for the light propagation from the
  // CLF to Celeste is approximately ~180 ns.

    fNsCorrection = 0;
    if (fCeleste == 1)
      fNsCorrection = 180;
    //--------------------------------------


    fCurEye = &eye; // acrobacy for TMinuit
    fCurEvent = &event; // acrobacy for TMinuit


    cout<<"Calling FindAxisHybrid (fHybridCore = 1)"<<"\n";
    //-----------------------------------------
    // This procedure will find the shower axis
    // minimizing the station-core distance.
    // (considering the station with the greater signal).
    // this step will help to avoid
    // finding a local minimum.
    fsignal_max = 90000000.0;
    fround = 1;
    fHybridCore = 1;
    if(this->FindAxisHybrid(eye))
        fFoundAxis=true;

    cout<<"Calling FindAxisHybrid (fHybridCore = 0)"<<"\n";



    //---------------------------------------------------
    // This procedure will find the shower axis
    // using the SD/FD times. The first guess
    // parameters are taking from the previous step.
    fHybridCore = 0;
    if(this->FindAxisHybrid(eye))
        fFoundAxis=true;

   

    // If the hybrid reconstruction failed try changing the
    // Chi0 first guess. Sometimes nearby inclined showers are
    // confused with upgoing showers. Reseting the first guess
    // helps correcting it.

    if ( (fChi0 <= 0. || fChi0 >  kPi ) && fStationAxisDistance > fDistanceStationSDP) {
      fChi0 = 90 * degree;
      fRp =  2000.;
      fHybridCore = 1;
      if(this->FindAxisHybrid(eye))
        fFoundAxis=true;

      fHybridCore = 0;
      if(this->FindAxisHybrid(eye))
        fFoundAxis=true;
    }



    //------------------------------------
    // Checking whether the station used for the hybrid
    // reconstruction is the right one to use


    cout<<"\n";
    cout<<"The total number of potential good stations is: "<<fTriggeredStations<<"\n";
    cout<<"(stations within "<<fDistanceStationSDP<<"m from the SDP)"<<"\n";
    cout<<"\n";

    //---------------------------------------
    //Checking whether the tank used for the
    //hybrid reconstruction gives a sencible
    //hybrid reconstruction.

    //The sd/fd offset (after the minization) has to
    //be smaller than 200ns.
    //The distance from the station to the SDP has
    //to be smaller than fDistanceStationSDPm.
    //The normilized time fit chisq has to be smaller than 50.

    // If any of these conditions is not satisfied, the program
    // will try the hybrid reconstruction with the next hottest
    // tank (if available).
    //----------------------------------------------------

    if (fTriggeredStations != 0) {
      cout<<"\n";
      cout<<"Tank with greatest signal:  "<<fHotStationId<<"\n";
      cout<<"   fSdHybridOffset "<<fSdHybridOffset<<"\n";
      cout<<"   fStationAxisDistance "<<fStationAxisDistance<<"\n";
      cout<<"   TimeFitChiSquare/dof "<<fChi2TimeFit<<" / "<<fNDofTimeFit<<"\n";
      cout<<"\n";
    } else {
      cout<<"\n";
      cout<<" There is not any station close to the \n";
      cout<<" SDP projection on the ground. Let's try \n";
      cout<<" another eye if available. Otherwise \n";
      cout<<" hybrid geometry finder has finished with this event \n";
      cout<<"\n";
      continue;
    }

    if (  fNDofTimeFit == 0 ) continue;  //This eye do not  have enough information
                                         // continue with other eye.


    bool GoodHybrid = true;
    if ( fabs(fSdHybridOffset) > 200. ||                //Cheking if the tested
         fStationAxisDistance > fDistanceStationSDP || //station makes a hybrid
         fChi2TimeFit/fNDofTimeFit > 50. )              //event
      GoodHybrid = false;
    else
      {
        nHybridStations++;
      }
    
         
    if ( nHybridStations == 0  && fround == fTriggeredStations  ) {
      cout<<" This event doesn't have any station consistent"<<" \n";
      cout<<" with the FD time and geometry."<<" \n";
      cout<<" from Eye: "<<eye.GetId()<<"\n";
      cout<<"\n";
      cout<<"No tank  ( last tank tested "<<fHotStationId
          <<") was consistent with the FD geometry and timing "<<"\n";
      cout<<"   fSdHybridOffset "<<fSdHybridOffset<<"\n";
      cout<<"   fStationAxisDistance "<<fStationAxisDistance<<"\n";
      cout<<"   TimeFitChiSquare/dof "<<fChi2TimeFit<<"/"<<fNDofTimeFit<<"\n";
      cout<<"   This was the trial number: "<<fround - 1<<"\n";
      continue;
    }    
    fsignal_max = fsignal_max_pre;
    
    //--------------------------------
    // Forcing to reconstruct with the first
    // and second hot station and choose  the
    // station that produces the smaller time
    // time fit chisqr.
    if ( fOptimize) {
      if ( GoodHybrid ) { 
        fChi0_test[0] = fChi0;
        fChi0Error_test[0] = fChi0Error;
        fT0_test[0] = fT0;
        fT0Error_test[0] = fT0Error;
        fRp_test[0] =fRp ;
        fRpError_test[0] = fRpError;
        fChi2TimeFit_test[0] = fChi2TimeFit;
        fNDofTimeFit_test[0] = fNDofTimeFit;
        fHotStationId_test[0] = fHotStationId;
        fround_test[0] = fround;
        for ( int i = 0; i < 3; i++)
          for ( int j = 0; j < 3; j++)
            fCovariance_test1[i][j] = fCovariance[i][j];
        
        fUseNthStation = 2;
      }
      else { 
        fChi2TimeFit_test[0] = -999.0;
        fUseNthStation = 1;
      }
    }
    
    //-------------------------------
    
    //-------------------------------------------------
    // the following if statement is to force the reconstruction
    // with the fUseNthStation hottest station (if any available)
    
    if (fround < fUseNthStation ) 
      GoodHybrid = false;
    //-----------------------------------------------

    if ( fUseNthStation > 1 &&  !fOptimize ) {
      cout<<"\n";
      cout<<"This hybrid reconstruction is forced to use the "<<fUseNthStation<<"the hottest station (if available)"<<"\n";
    } 

    if ( fUseNthStation > 1 &&  fOptimize ) {
      cout<<"\n";
      cout<<"Now checking the second hot station (if available)"<<"\n";
    } 

       
    //Cheking if other stations can make a hybrid event
    while ( !GoodHybrid && fTriggeredStations > fround )
      {
        fround++;
        if( fDebuging == 1) {
          cout << "\n";
          cout << "The tank with greatest signal ("<<fHotStationId<<") "<<"\n";
          cout << "the hybrid reconstruction: "<<"\n";
          cout << "   fSdHybridOffset "<<fSdHybridOffset<<"\n";
          cout << "   fStationAxisDistance "<<fStationAxisDistance<<"\n";
          cout << "   TimeFitChiSquare/dof "<<fChi2TimeFit<<" / "<<fNDofTimeFit<<"\n";
          cout << "We will try the "<<"\n";
          cout << "reconstruction with the next tank with the greatest \n";
          cout << "signal. This is the trial number: "<<fround<<"\n";
        }
        //-----------------------------------------
        // This procedure will find the shower axis
        // minimizing the station-core distance.
        // (considering the station with the greater signal).
        // this step will help to avoid
        // finding a local minimum.
        fHybridCore = 1;
        cout<<"Calling FindAxisHybrid (fHybridCore = 1)"<<" round:" <<fround<<"\n";
        // fround tells FindAxisHybrid whether to use the hottest
        // or the second , or the third (etc) hottest tank.
        if(this->FindAxisHybrid(eye))
            fFoundAxis=true;


        //---------------------------------------------------
        // This procedure will find the shower axis
        // using the SD/FD times. The first guess
        // parameters are taking from the previous step.

        fHybridCore = 0;
        cout<<"Calling FindAxisHybrid (fHybridCore = 0)"<<" round:"<<fround<<"\n";
        if(this->FindAxisHybrid(eye))
          fFoundAxis=true;


        // If the hybrid reconstruction failed try changing the
        // Chi0 first guess. Sometimes near by inclined showers are
        // confused with upgoing showers. Reseting the first guess
        // helps correcting it.

        if ( (fChi0 <= 0. || fChi0 >  kPi )  && fStationAxisDistance > fDistanceStationSDP) {
          fChi0 = 90 * degree;
          fRp =  2000.;
          fHybridCore = 1;
          if(this->FindAxisHybrid(eye))
            fFoundAxis=true;
          fHybridCore = 0;
          if(this->FindAxisHybrid(eye))
            fFoundAxis=true;
        }
        
        //----------------------------
        
        if ( fabs(fSdHybridOffset) > 200 ||
             fStationAxisDistance > fDistanceStationSDP ||
             fChi2TimeFit/fNDofTimeFit > 50 ) {
          cout<<" Station "<<fHotStationId<<" is not consistent \n";
          cout<<" with the FD event \n";
          fsignal_max = fsignal_max_pre;
        }
        else  {
          cout<<"The tank with greatest signal ("<<fHotStationId<<") "<<"\n";
          cout<<"the hybrid reconstruction: "<<"\n";
          cout<<"   fSdHybridOffset "<<fSdHybridOffset<<"\n";
          cout<<"   fStationAxisDistance "<<fStationAxisDistance<<"\n";
          cout<<"   TimeFitChiSquare/dof "<<fChi2TimeFit<<" / "<<fNDofTimeFit<<"\n";
          cout<<"\n";
          nHybridStations++;
          if (fround < fUseNthStation ) {
            cout<<" this is a good station ("<<fHotStationId <<"), but we will force the  \n";
            cout<<" reconstruction using the "<<fUseNthStation<<"th hottest station. \n";
            
            fsignal_max = fsignal_max_pre;
            continue;
          }
          else {
            GoodHybrid = true;
            break;
          }
        }
        if ( !GoodHybrid  && fround == fTriggeredStations  ) {
          cout<<" This event doesn't have any station consistent"<<" \n";
          cout<<"  with the FD time and geometry."<<" \n";
          cout<<"\n";
          cout<<"No tank  ( last tank tested "<<fHotStationId
              <<") was consistent with the FD geometry and timing "<<"\n";
          cout<<"   fSdHybridOffset "<<fSdHybridOffset<<"\n";
          cout<<"   fStationAxisDistance "<<fStationAxisDistance<<"\n";
          cout<<"   TimeFitChiSquare/dof "<<fChi2TimeFit<<"/"<<fNDofTimeFit<<"\n";
          cout<<"This was the trial number: "<<fround - 1<<"\n";
          break;
        }
      }


    //=====================================
   if ( fOptimize) {
      
      if ( GoodHybrid ) { 
        fChi0_test[1]        = fChi0;
        fChi0Error_test[1]   = fChi0Error;
        fT0_test[1]          = fT0;
        fT0Error_test[1]     = fT0Error;
        fRp_test[1]          = fRp;
        fRpError_test[1]     = fRpError;
        fChi2TimeFit_test[1] = fChi2TimeFit;
        fNDofTimeFit_test[1] = fNDofTimeFit;
        fHotStationId_test[1] = fHotStationId;
        fround_test[1] = fround;
        for ( int i = 0; i < 3; i++)
          for ( int j = 0; j < 3; j++)
            fCovariance_test2[i][j] = fCovariance[i][j];

      }
      else { 
        fChi2TimeFit_test[1] = -999.0;
      }
    }
    
    //------------------------------------
    // Selecting which station to use
    if ( nHybridStations == 0 )  continue;

    if (nHybridStations == 1) 
      if ( fChi2TimeFit_test[0] > 0. ) FirstGoodSt = true;

    if (nHybridStations == 2) {
      if ( ( fChi2TimeFit_test[0]/fNDofTimeFit_test[0] <  
             (fChi2TimeFit_test[1]/fNDofTimeFit_test[1]) + 1.0 ) ) 
        FirstGoodSt = true;         
    }
    

    if (FirstGoodSt) 
      UseStation = 0;
    else
      UseStation = 1;
   

    //-------------------------------------------------------------------------------
    if ( fOptimize ) {

      fChi0         = fChi0_test[UseStation];
      fChi0Error    = fChi0Error_test[UseStation];
      fT0           = fT0_test[UseStation];
      fT0Error      = fT0Error_test[UseStation];
      fRp           = fRp_test[UseStation];
      fRpError      = fRpError_test[UseStation];
      fChi2TimeFit  = fChi2TimeFit_test[UseStation];
      fNDofTimeFit  = fNDofTimeFit_test[UseStation];
      fHotStationId  = fHotStationId_test[UseStation];
      fround         = fround_test[UseStation];
      for ( int i = 0; i < 3; i++)
        {
          for ( int j = 0; j < 3; j++)
            {
              if ( UseStation == 0 )
                fCovariance[i][j] = fCovariance_test1[i][j];
              else
                fCovariance[i][j] = fCovariance_test2[i][j];
            }
        }
      
      
      
      if ( nHybridStations == 2 &&
           fabs( fChi0_test[1]-fChi0_test[0] ) > 1.0 * degree &&
           fabs(  fChi2TimeFit_test[0] - fChi2TimeFit_test[1] ) < 0.5 ) {
        cout<<" There are two ambigous geometries (similar time fit chisqr) \n";
        cout<<" with a difference greater than 1.0 deg. \n";
        cout<<" The profile reconstruction may help two identify \n";
        cout<<" the correct geometry. This event  should be \n";
        cout<<" FLAGGED for further quality check. \n";
      }
      
    }
    
    
    //---------------------------------------
    // Readmitting pixels according to the hybrid geometry
    // and then rejecting outliers one by one
    while( ReadmitPixel(eye) )
      if(this->FindAxisHybrid(eye))  fFoundAxis=true;
    

    while ( RejectPixel(eye) )
      if(this->FindAxisHybrid(eye))  fFoundAxis=true;
    //-----------------------------------------------------

    
    if(!fFoundAxis) {
      cerr << " FATAL ERROR in HybridGeometryFinder::FindAxisHybrid  Minuit MIGRAD failed (Eye: "<<eye.GetId()<<"), try other Eye if available"<<endl;
      continue;
    }
    
    cout<<"  number of pixels remaining "<<eyerecdata.GetNTimeFitPixels()<<"\n";    cout<<"\n";

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

    // catch some very rare crashs
    if (fChi0 == 0 && fRp == 0 && fT0 == 0)
      continue;

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

    eyerecdata.SetChiZero(fChi0,fChi0Error);
    eyerecdata.SetTZero(fT0,fT0Error);
    eyerecdata.SetRp(fRp,fRpError);
    eyerecdata.SetTimeFitChiSquare(fChi2TimeFit, fNDofTimeFit);

    cout<<"Timing-Chisqr / dof "<<fChi2TimeFit<<" / "<<fNDofTimeFit<<"\n";
    cout<<"\n";
    cout << " Chi0=" << fChi0/deg << " (" << fChi0Error/deg << ") deg "
         << " Rp=" << fRp/m << " (" << fRpError/m << ") m "
         << " T0=" << fT0/ns << " (" << fT0Error/ns << ") ns ";
    cout<<"\n";

    ostringstream info;
    info << " Chi0=" << fChi0/deg << " (" << fChi0Error/deg << ") deg "
         << " Rp=" << fRp/m << " (" << fRpError/m << ") m "
         << " T0=" << fT0/ns << " (" << fT0Error/ns << ") ns ";
    INFO(info);
    // correlation coefficiens

    const double rRpT0   = fCovariance[2][1]/fRpError/fT0Error;
    const double rRpChi0 = fCovariance[2][0]/fRpError/fChi0Error;
    const 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(fChi2TimeFit, fNDofTimeFit);
      telRecData.SetTimeFitCorrelations(rRpT0, rRpChi0, rChi0T0);
    }

    //==========================================
    // calculate from chi0, Rp the axis parameters.
    Vector sdp = eyerecdata.GetSDP();
    Vector vertical2(0,0,1,eyeCS);
    Vector horizontalInSDP = cross(sdp,vertical2); // horizontal
    horizontalInSDP.Normalize();
    Transformation rot (Transformation::Rotation(-fChi0,sdp, eyeCS));
    Vector axis (rot* horizontalInSDP);             // apply rotation
    axis.Normalize();
    CoordinateSystemPtr siteCS =
      det::Detector::GetInstance().GetSiteCoordinateSystem();
    Vector core_eye_vec = fRp / sin( kPi - fChi0) * horizontalInSDP;
    const Point& eyeposition =
      det::Detector::GetInstance().GetFDetector().GetEye(eye).GetPosition();
    Point core = eyeposition + core_eye_vec;
    //Setting the Axis parameters
    if (!event.HasRecShower()) event.MakeRecShower();
    ShowerRecData& recshower = event.GetRecShower();
    recshower.SetAxis(axis);
    recshower.SetCorePosition(core);
    if (!eyerecdata.HasFRecShower())
      eyerecdata.MakeFRecShower();
    eyerecdata.GetFRecShower().SetAxis(axis);
    eyerecdata.GetFRecShower().SetCorePosition(core);
    eyerecdata.GetFRecShower().SetCoreTime(fCurEye->GetHeader().GetTimeStamp()-
                                           TimeInterval(1000*100*ns),
                                           fRp, fChi0, fT0);
    eyerecdata.GetFRecShower().SetSDTimeResidual(fSdHybridOffset);
    eyerecdata.GetFRecShower().AddStationId(fHotStationId);

  }// for eyeiter


  return eSuccess;
}


bool
HybridGeometryFinder::FindAxisHybrid(fevt::Eye& /*eye*/)
{
  const int npar = 3; // number of parameters

   TMinuit theMinuit(npar); // 3 parameter fit

   theMinuit.SetPrintLevel(-1); // minuit verbosity

   theMinuit.SetFCN(HybridGeometryFinder::MinuitFitFuncHybrid);

   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



   vstart[0] = fChi0 ;
   vstart[1] = fT0;
   vstart[2] = fRp ;

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


   //init function
   arglist[0] = 1;
   theMinuit.mnexcm("CALI", arglist, 1, ierflag);

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

   if (ierflag) {
       ERROR ("Minuit MIGRAD failed");
       return false;
   }

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

   return true;
}


void
HybridGeometryFinder::MinuitFitFuncHybrid(int& /*npar*/, double* const /*gin*/,
                                          double& value, double* const par,
                                          const int iflag)
{
  fevt::EyeRecData& recdata = fCurEye->GetRecData();

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

  double sum_charge = 0;
  double chisq_time = 0;
  double chisq_hybrid = 0;

  const fdet::FTimeFitModel& timeFitModel = fdet::FTimeFitModel::GetInstance();

  static vector<double> t_i;
  static vector<double> t_i_Err;
  static vector<double> chi_i;
  static vector<int> telID;
  //static vector<double> charge;   not used
  static int eyeID;
  static CoordinateSystemPtr eyeCS;
  static Point eyeposition;
  static Vector horizontalWithinSDP;

  if (iflag == 1) {  // fnc init
    t_i.clear();
    t_i_Err.clear();
    chi_i.clear();
    telID.clear();

    eyeID = fCurEye->GetId();
    eyeCS = det::Detector::GetInstance().GetFDetector().GetEye(*fCurEye).GetEyeCoordinateSystem();

    //get pixel data
    for (EyeRecData::ConstPixelIterator pixelIter= recdata.TimeFitPixelsBegin();
         pixelIter!=recdata.TimeFitPixelsEnd(); ++pixelIter){
      const PixelRecData& pixrec = pixelIter->GetRecData();
      telID.push_back(pixelIter->GetTelescopeId());

      t_i.push_back(pixrec.GetT_i().GetNanoSecond());
      t_i_Err.push_back(pixrec.GetT_iError().GetInterval());
      chi_i.push_back( pixrec.GetChi_i());
      // charge.push_back(pixrec.GetTotalCharge());


    }

    const Vector vertical(0,0,1,eyeCS,Vector::kCartesian);
    horizontalWithinSDP = cross(fSDP, vertical);
    horizontalWithinSDP /= horizontalWithinSDP.GetMag();

    eyeposition=det::Detector::GetInstance().GetFDetector().GetEye(*fCurEye).GetPosition();
  }

  const double sdpTheta = fSDP.GetTheta(eyeCS);
  for(unsigned int pix=0; pix<t_i.size(); ++pix){
    const double t_expected = timeFitModel.GetTimeAtAperture(t_0, r_p, chi_0, chi_i[pix], sdpTheta, eyeID, telID[pix]);
    const double tmp = t_i[pix] - t_expected;
    //  NOT weighting with the pixel charge
    const double chargeOne = 1.0;
    chisq_time += tmp*tmp *chargeOne / (t_i_Err[pix] * t_i_Err[pix]);
    sum_charge +=chargeOne;
  }
  chisq_time /= sum_charge;

  static vector<Point>  stationpos;
  static vector<int>    stationId;
  static vector<double> station_d_to_SDP;
  static vector<double> stNanosecond;
  static vector<double> stSignal;

  if (iflag == 1) {  // fnc init

    stationpos.clear();
    stationId.clear();
    station_d_to_SDP.clear();
    stNanosecond.clear();
    stSignal.clear();

    sevt::SEvent& sevent = fCurEvent->GetSEvent();
    for(SEvent::StationIterator sditer= sevent.StationsBegin();
        sditer != sevent.StationsEnd(); ++sditer){
      if (!sditer->HasRecData())
        continue;
      sevt::Station& station= *sditer;
      if (!station.IsCandidate() && !station.IsSilent() )
        continue;
      const sdet::Station& detstation =        det::Detector::GetInstance().GetSDetector().GetStation(station);
      if(detstation.IsDense())
        continue;

      //Estimating distance to the SDP
      const Vector eye_station = detstation.GetPosition() - eyeposition;
      double angle = Angle(eye_station,fSDP);
      double distToSDP= abs(eye_station.GetMag() * cos(angle));

      double nanosecond = station.GetRecData().GetSignalStartTime().GetGPSNanoSecond();
      // Checking if the nanoseconds correspond to the same second in
      // in the SD and FD times
      int SDSecond = station.GetRecData().GetSignalStartTime().GetGPSSecond();
      int FDSecond = fCurEye->GetHeader().GetTimeStamp().GetGPSSecond();
      if ( SDSecond > FDSecond ) nanosecond += 1000000000;
      if ( SDSecond < FDSecond ) nanosecond -= 1000000000;
      double signal = station.GetRecData().GetTotalSignal();

      stationpos.push_back(detstation.GetPosition());
      stationId.push_back(station.GetId());
      station_d_to_SDP.push_back(distToSDP);
      stNanosecond.push_back(nanosecond);
      stSignal.push_back(signal);

    }
  }

  double signal_max = 0;
  int hotStationIter = 0;
  int NofStations = 0;
  fTriggeredStations = 0;
  
  for(unsigned int statIt=0; statIt<stationpos.size(); ++statIt){

    if (station_d_to_SDP[statIt] < fDistanceStationSDP)
      fTriggeredStations++;
    else
      continue;
    
    if (fround == 1) {
      if (stSignal[statIt] > signal_max ) {
        signal_max = stSignal[statIt];
        hotStationIter = statIt;
        NofStations++;
        fsignal_max = stSignal[statIt];
        fsignal_max_pre = signal_max;
        fHotStationId = stationId[statIt];
      }
    }
    else {
      if (stSignal[statIt] > signal_max && stSignal[statIt] < fsignal_max){
        signal_max = stSignal[statIt];
        fsignal_max_pre = signal_max;
        hotStationIter = statIt;
        NofStations++;
        fHotStationId = stationId[statIt];
      }
    }
  }//for stations


  //get shower axis
  Transformation rot(Transformation::Rotation(-chi_0, fSDP, eyeCS));
  Vector axis(rot* horizontalWithinSDP);
  axis /= axis.GetMag();


  const double nx = axis.GetX(eyeCS);
  const double ny = axis.GetY(eyeCS);
  const double nz = axis.GetZ(eyeCS);

  const Vector core =  horizontalWithinSDP*( r_p/ sin(kPi-chi_0));
  const double tcore =  t_0 + r_p/(kSpeedOfLight * tan(kPi - chi_0));
  const Point PointCore = eyeposition + core;

  int stations_used = 0;
  if (NofStations > 0 ) {
    const  double x_hot = stationpos[hotStationIter].GetX(eyeCS);

    //---------------------------------
    UTMPoint UTMposHot(stationpos[hotStationIter], ReferenceEllipsoid::GetWGS84());

    double sum_total = 0;
    for(unsigned int stiter=0; stiter<stationId.size(); ++stiter){

      //--------------------------------------
      //  This statement is to consider the
      //  stations around the hottest station
      //  in the hybrid chisq.

      // 500 means that only the hottest station
      // will be used for the hybrid reconstruction.
      const Vector distance_to_hot= stationpos[stiter]- stationpos[hotStationIter];

      if ( distance_to_hot.GetMag() > 500 )
        continue;

      const Vector st_core= stationpos[stiter] - PointCore;
      const Vector sh_axis ( nx, ny, nz, eyeCS, Vector::kCartesian);
      const double sh_axis_VS_st_core_angle = Angle(st_core,sh_axis);
      const double delta_d = st_core.GetMag() * std::cos(sh_axis_VS_st_core_angle);

      //Estimating station distance to the Axis
      fStationAxisDistance = st_core.GetMag() * std::sin(sh_axis_VS_st_core_angle);

      double station_to_axis = st_core.GetMag() * std::sin(sh_axis_VS_st_core_angle);
      double exp_t_st = tcore - delta_d/kSpeedOfLight;

      // time delay due to the shower front curvature
      const double d2 = station_to_axis * station_to_axis;
      const double curvature = fVerticalCurvature*sh_axis.GetCosTheta(eyeCS);
      const double timeDelay = curvature * d2 / (2*kSpeedOfLight);
      exp_t_st += timeDelay;

      // Get event TimeStamp

      utl::TimeStamp time2 = fCurEye->GetHeader().GetTimeStamp();
      // double(time2.GetGPSNanoSecond()) gives the time at the end of the FD FADC trace.
      // by substracting 1000 *100 ns we get the start time of  the FD FADC trace.
      // exp_t_st, is the expected trigger time of the station.
      // This time is currently relative to the start of the FD FADC trace. We'll add 'time2.GetGPSNanoSecond()'
      // (the eye trigger  time' to get the absolute time and compare with the observed station time.
      double eye_trigger_time = double(time2.GetGPSNanoSecond());

      //-------------------------
      eye_trigger_time -= 1000 *100;
      exp_t_st += eye_trigger_time;

      double SD_Time = stNanosecond[stiter]; //nanosecond;
      double FD_Time = exp_t_st + fNsCorrection;
      double chi_hybrid = SD_Time - FD_Time;
      double hybrid_error = 100;

      if (fHybridCore == 1) {
        //Minimizing the hot_station-core distance (first step)
        const double x_st=stationpos[stiter].GetX(eyeCS);
        if (abs(x_hot -x_st) < 1 ) {
          chisq_hybrid +=  st_core.GetMag()* st_core.GetMag()/(10*10);
          sum_total++;
          stations_used++;
          break;
        }
      }
      else {
        //------------------------------------------------
        //Minimizind the SD/FD time offset (second step).
        // A different weight is given to different stations.
        // The weight is proportional to the inverse of the
        // distance from the station to the core.

        const double hybrid_weight = (1/station_to_axis) * (1/station_to_axis);
        chisq_hybrid += chi_hybrid * chi_hybrid * hybrid_weight / (hybrid_error*hybrid_error);
        sum_total += hybrid_weight;
        stations_used++;
        fSdHybridOffset = chi_hybrid;
      }
    }
    chisq_hybrid /= sum_total;
  }
  else {
    chisq_hybrid = 0;
  }

  if( fDebuging == 1)
    cout<<"chisq time: "<<chisq_time<<" chisq hybrid: "<<chisq_hybrid<<" \n";
  //================================= SD (vitor) (end)

  chisq = 0;
  // Normalized hybrid chisq times the number of stations.
  chisq_hybrid *= stations_used;
  // Normalized timing chisq times the number of pixels.
  chisq_time *= recdata.GetNTimeFitPixels();

  //Total chisq.
  chisq += chisq_time;

  chisq += chisq_hybrid;
  if( fDebuging == 1)
    cout<<"chisq time (TOTAL): "
        <<chisq_time<<" chisq hybrid: "
        <<chisq_hybrid<<" total: "
        <<chisq<<" \n";

  // Minuit needs the Total (not the reduced) chisq to estimate the errors
  // correctly
  value = chisq;

  fChi2TimeFit = chisq_time;
  fNDofTimeFit = recdata.GetNTimeFitPixels()-3;

}


bool
HybridGeometryFinder::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++;
      cout << "readmitting pixel  " << evtpixel.GetId() << " for the Timing Fit\n";
    }

  }
  return ReadmittedPixels > 0;
}


bool
HybridGeometryFinder::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;
    }
  }

  if (RejectedPixels>0)
    return true;
  else
    return false;
}


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


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