G/PCGFitter.cc

Go to the documentation of this file.

#define DEBUG 1
#include "PCGFitter.h"

#include <utl/Point.h>
#include <utl/Vector.h>
#include <utl/CoordinateSystemPtr.h>
#include <utl/Transformation.h>
#include <utl/TimeInterval.h>
#include <utl/AugerUnits.h>
#include <utl/PhysicalConstants.h>
#include <utl/TabulatedFunctionErrors.h>

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

#include <fevt/Telescope.h>
#include <fevt/TelescopeRecData.h>

#include <det/Detector.h>
#include <fdet/FDetector.h>
#include <fdet/Eye.h>
#include <fdet/Pixel.h>

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

#include "ChiZeroRegression.h"
#include "ProfileChi2.h"

#include "../../General/RecDataWriterNG/ConversionUtil.h"

#include <limits>

#define PCGF_DEBUGLOG(x) { cerr << "DEBUGLOG " << x << endl;}


using namespace std;
using namespace fevt;
using namespace utl;

namespace FdProfileConstrainedGeometryFitPG {

void
PCGFitter::Init()
{
  Branch topBranch = fwk::CentralConfig::GetInstance()->GetTopBranch("FdProfileConstrainedGeometryFitPG");

  string apLightOption;
  topBranch.GetChild("apLightMethod").GetData(apLightOption);
  fApLightMethod = (apLightOption=="eExternal") ?
    eExternal :
    (apLightOption=="eInternal")?
    eInternal : eCallForEach;

  topBranch.GetChild("checkUnderground").GetData(fCheckUnderground);
  topBranch.GetChild("prescan").GetData(fPrescan);
  topBranch.GetChild("scanOnly").GetData(fScanOnly);
  topBranch.GetChild("scanStep").GetData(fScanStep);
  topBranch.GetChild("scanStart").GetData(fScanStart);
  topBranch.GetChild("scanStop").GetData(fScanStop);
  topBranch.GetChild("useLightFlux").GetData(fUseLightFlux);
  topBranch.GetChild("skipNegativeT0").GetData(fSkipNegativeT0);
  topBranch.GetChild("delZeroLightFlux").GetData(fDelZeroFlux);
  topBranch.GetChild("TimeFitModel").GetData(fRealAtm);
  topBranch.GetChild("TimeFitDeexcitation").GetData(fDeex);
  topBranch.GetChild("emissionPointCorrection").GetData(fEmissionPointCorrection);
  ostringstream info;
  info << "\n"
          "\t apLightMethod is " << apLightOption << "\n"
          "\t checkUnderground is " << (fCheckUnderground?"on":"off") << "\n"
          "\t prescan is " << (fPrescan?"on":"off") << "\n"
          "\t scanOnly is " << (fScanOnly?"on":"off") << "\n"
          "\t scanStep: " << fScanStep/degree << " deg\n"
          "\t scanStart: " << fScanStart/degree << " deg\n"
          "\t scanStop: " << fScanStop/degree << " deg\n"
          "\t delZeroLightFlux is " << (fDelZeroFlux?"on":"off") << "\n"
          "\t useLightFlux is " << (fUseLightFlux?"on":"off") << "\n"
          "\t skipNegativeT0 is " << (fSkipNegativeT0?"on":"off") << "\n"
          "\t emissionPointCorrection is " << (fEmissionPointCorrection?"on":"off") << "\n"
          "\t TimeFitModel is " << (fRealAtm?"RealisticAtm":"VacuumAtm") << "\n"
          "\t TimeFitDeexcitation is " << (fDeex?"on":"off") << "\n";
  INFO(info);

  fPrecise = false;
  if (fApLightMethod == eCallForEach) {
    fApertureLightFinder.Init();
    fPrecise = true;
  } else if (fApLightMethod == eInternal) {
    fApLight.Init();
  }

  fProfileFit.Init();

  fProfileFit.SetUseLightFlux(fUseLightFlux);
}


inline
bool
PCGFitter::Underground(const double chi0, const double rp, const Eye &eye)
{
  if (!fCheckUnderground)
    return false;

  const EyeRecData& eyeRec = eye.GetRecData();

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

  const Point telPos = detFD.GetEye(eye.GetId()).GetTelescope(fTelId).GetPosition();
  utl::CoordinateSystemPtr telCS = fwk::LocalCoordinateSystem::Create(telPos);
  const Vector sdp = eye.GetTelescope(fTelId).GetRecData().GetSDP();
  const Vector vertical(0,0,1, telCS);
  const Vector horizontalInSDP = Normalized(Cross(sdp, vertical));

  const Transformation rot = Transformation::Rotation(-chi0, sdp, telCS);
  Vector axis = Normalized(rot*horizontalInSDP);

  const double coreDistance = rp / sin(chi0);
  const Vector coreTelVec = coreDistance * horizontalInSDP;
  const Point core = telPos + coreTelVec;

  for (EyeRecData::ConstPixelIterator pix = eyeRec.TimeFitPixelsBegin(),
       end = eyeRec.TimeFitPixelsEnd(); pix != end; ++pix) {

    int telId = pix->GetTelescopeId();
    //if (telId == fTelId) {

      const fdet::Telescope& curTel = detFD.GetEye(eye.GetId()).GetTelescope(telId);
      const Point curtelPos = curTel.GetPosition();
      const CoordinateSystemPtr curtelCS = fwk::LocalCoordinateSystem::Create(curtelPos);

      const Point core_in_curtel = core-axis*core.GetZ(curtelCS)/cos(axis.GetTheta(curtelCS));

      const Vector curtel_to_core = core_in_curtel-curtelPos;
      double chi0_in_curtel = Angle(curtel_to_core, axis);
      if (axis.GetTheta(curtelCS) > kPi/2.) {
        chi0_in_curtel = -(kPi-chi0_in_curtel);
      }
      const utl::AxialVector SDP_in_curtel = Normalized(Cross(axis, curtel_to_core));

      const Vector vertical_in_curtel(0,0,1, curtelCS);
      const Vector horizontalInSDP_in_curtel = Normalized(Cross(SDP_in_curtel, vertical_in_curtel));

      const Vector pixelDirection = curTel.GetPixel(pix->GetId()).GetDirection();
      const Vector pixelDirInSDP = pixelDirection - (SDP_in_curtel * pixelDirection)*SDP_in_curtel;
      const double chii = Angle(pixelDirInSDP, horizontalInSDP_in_curtel);
      //following is too simple! - eyeCopy does not make new fevt::Pixel and in RecData only pointers are stored, this makes mess
      //const PixelRecData& pixRec = pix->GetRecData();
      //const double chii = pixRec.GetChi_i();

      //cout << "chii " << chii << " pixel " << (*pix).GetId() << " tel "<<telId<<" eye "<<eye.GetId()<<"\n";
      if (sin(chi0_in_curtel - chii) < 0) {
        PCGF_DEBUGLOG("Under ground");
        return true;
      }
    //}
  }
  return false;
}


double
PCGFitter::operator()(const vector<double>& par)
{
  const double chi0 = par[0];

  double rp = 0;
  double t0 = 0;
  double rpErr = 0;
  double t0Err = 0;
  const double chi2 = CombinedChi2(chi0, rp, rpErr, t0, t0Err, PCGFData::eFit, true); 
  return chi2; 
}


double
PCGFitter::CombinedChi2(const double chi0, double &rp, double &rpErr, double &t0, double &t0Err,
                        const fevt::PCGFData::EPCGFStatus status,
                        const bool verbose)
{
  const Eye& eye = (*fEye);

  //horizontal shower with undefined core
  if (fabs(chi0-180*deg) < 1e-4*deg || fabs(chi0+180*deg) < 1e-4*deg || fabs(chi0-0*deg) < 1e-4*deg || fabs(chi0-360*deg) < 1e-4*deg)
    return numeric_limits<double>::infinity();

  double chi2timing; // rp, t0;
  const ChiZeroRegression& chi0reg = *fChi0Regression;
  chi0reg(chi0, chi2timing, rp, t0);

  if (fSkipNegativeT0 && (t0 < 0))
    return numeric_limits<double>::infinity();

  if (verbose)
    cout << "Chi0: " << chi0/degree << "\t"
            "Rp: " << rp/m << "\t"
            "T0: " << t0/ns << "\t "
            "Chi^2_t: " << chi2timing << "\n";

  //Check after regression: fill rp, t0 correctly anyway; probably never needed ...
  if (Underground(chi0, rp, eye))
    return numeric_limits<double>::infinity();

  evt::Event eventCopy;
  Eye& eyeCopy = CopyEye(eye, eventCopy);

  FillParams(eyeCopy, chi0, rp, t0);
  if (!AdjustGeometry(eyeCopy)) {
    //ERROR("could not adjust geometry");
    return numeric_limits<double>::infinity();
  }
  
  if (fApLightMethod == eCallForEach) {
    const bool haveLight =
      fApertureLightFinder.Run(eventCopy) == fwk::VModule::eSuccess && eyeCopy.GetRecData().HasLightFlux();
    if (!haveLight)
      return numeric_limits<double>::infinity();

    if (fDelZeroFlux) {
      //VN remove all light flux bins comatible with zero
      TabulatedFunctionErrors& photonTrace = eyeCopy.GetRecData().GetLightFlux(fevt::FdConstants::eTotal);
      TabulatedFunctionErrors& bgTrace = eyeCopy.GetRecData().GetLightFlux(fevt::FdConstants::eBackground);
      vector<double> phX, phXerr, phY, phYerr;
      vector<double> bgX, bgXerr, bgY, bgYerr;
      for (unsigned int i = 0; i < photonTrace.GetNPoints(); ++i) {
        if ((photonTrace.GetY(i)-photonTrace.GetYErr(i)) > 0.) {
          phX.push_back(photonTrace.GetX(i));
          phXerr.push_back(photonTrace.GetXErr(i));
          phY.push_back(photonTrace.GetY(i));
          phYerr.push_back(photonTrace.GetYErr(i));
          bgX.push_back(bgTrace.GetX(i));
          bgXerr.push_back(bgTrace.GetXErr(i));
          bgY.push_back(bgTrace.GetY(i));
          bgYerr.push_back(bgTrace.GetYErr(i));
        }
      }
      photonTrace.Clear();
      bgTrace.Clear();
      for (unsigned int i = 0; i < phX.size(); ++i) {
        photonTrace.PushBack(phX[i], phXerr[i], phY[i], phYerr[i]);
        bgTrace.PushBack(bgX[i], bgXerr[i], bgY[i], bgYerr[i]);
      }

      for (fevt::Eye::TelescopeIterator telIter =
           eyeCopy.TelescopesBegin(fevt::ComponentSelector::eInDAQ);
         telIter != eyeCopy.TelescopesEnd(fevt::ComponentSelector::eInDAQ);
         ++telIter) {
        if (!telIter->HasRecData()) continue;
        TelescopeRecData& telRecData = telIter->GetRecData();
        if (!telRecData.HasLightFlux()) continue;
        TabulatedFunctionErrors& photonTrace = telRecData.GetLightFlux(fevt::FdConstants::eTotal);
        TabulatedFunctionErrors& bgTrace = telRecData.GetLightFlux(fevt::FdConstants::eBackground);
        vector<double> phX, phXerr, phY, phYerr;
        vector<double> bgX, bgXerr, bgY, bgYerr;
        for (unsigned int i = 0; i < photonTrace.GetNPoints(); ++i) {
          if (photonTrace.GetY(i)-photonTrace.GetYErr(i) > 0.) {
            phX.push_back(photonTrace.GetX(i));
            phXerr.push_back(photonTrace.GetXErr(i));
            phY.push_back(photonTrace.GetY(i));
            phYerr.push_back(photonTrace.GetYErr(i));
            bgX.push_back(bgTrace.GetX(i));
            bgXerr.push_back(bgTrace.GetXErr(i));
            bgY.push_back(bgTrace.GetY(i));
            bgYerr.push_back(bgTrace.GetYErr(i));
          }
        }
        photonTrace.Clear();
        bgTrace.Clear();
        for (unsigned int i = 0; i < phX.size(); ++i) {
          photonTrace.PushBack(phX[i], phXerr[i], phY[i], phYerr[i]);
          bgTrace.PushBack(bgX[i], bgXerr[i], bgY[i], bgYerr[i]);
        }
      }
    }
  }

  evt::Event eventCopy2;
  Eye& eyeCopy2 = CopyEye(eyeCopy, eventCopy2);
  if (std::isfinite(fProfileFit(eyeCopy2, fPrecise))) {
    double Xmax = eyeCopy2.GetRecData().GetFRecShower().GetGHParameters().GetXMax();
    chi0reg(chi0, chi2timing, rp, rpErr, t0, t0Err, Xmax);

    FillParams(eyeCopy, chi0, rp, t0);
    if (!AdjustGeometry(eyeCopy)) {
      return numeric_limits<double>::infinity();
    }
  }

  const double chi2profile = fProfileFit(eyeCopy, fPrecise);
  const double constraint1 = fProfileFit.GetShapeConstraint1Chi2();
  const double constraint2 = fProfileFit.GetShapeConstraint2Chi2();
  const double constraint3 = fProfileFit.GetShapeConstraint3Chi2();
  if (verbose) {
    cout << "Chi^2_prof: " << chi2profile << ",  including constraints 1/2/3: " << constraint1 << "/" << constraint2 << "/" <<  constraint3 << '\n';
    if (fEmissionPointCorrection)
      cout << "Chi^2_time after emission point correction: " << chi2timing << '\n';
  }
  
  map<string, double> mapchi2;
  mapchi2["time"] = chi2timing;
  mapchi2["prof"] = chi2profile-constraint1-constraint2-constraint3;//sqrt(pow(chi2profile,2) - pow(constraint1,2) - pow(constraint2,2) - pow(constraint3,2) );
  mapchi2["constraint1"] = constraint1;
  mapchi2["constraint2"] = constraint2;
  mapchi2["constraint3"] = constraint3;
  map<string, double> mapNDof;
  mapNDof["time"] = chi0reg.GetNDof();
  mapNDof["prof"] = fProfileFit.GetNDof();
  mapNDof["constraint1"] = 1;
  mapNDof["constraint2"] = 1;
  mapNDof["constraint3"] = 1;
  fPCGFData.push_back(PCGFData(chi0, rp, t0, status, mapchi2, mapNDof));
  
  return chi2profile+chi2timing;
}


bool
PCGFitter::Run(Eye &eye)
{
  fEye = &eye;

  //VN to manage static fitter
  fProfileFit.SetFitterPars();

  fPCGFData.clear();

  // keep initial value
  {
    const EyeRecData& eyeRecData = eye.GetRecData();
    const double t0 = eyeRecData.GetTZero();
    const double chi0 = eyeRecData.GetChiZero();
    const double rp = eyeRecData.GetRp();
    
    map<string, double> mapchi2;
    mapchi2["time"] = eyeRecData.GetTimeFitChiSquare();
    map<string, double> mapNDof;
    mapNDof["time"] = eyeRecData.GetTimeFitNDof();
    // typically we have no profile fit at this stage...
    /*
    ShowerFRecData& recShower = eyeRecData.GetFRecShower();
    mapchi2["prof"] = sqrt(recShower.GetGHParameters().GetChiSquare();
    //recShower.GetGHParameters().GetNdof()) :
    */
    fPCGFData.push_back(PCGFData(chi0, rp, t0, PCGFData::ePreFit, mapchi2, mapNDof));
  }

  ChiZeroRegression chi0reg(eye, fTelId);
  if (fTelId == 0) {
    return false;
  }
  chi0reg.SetUseLightFlux(fUseLightFlux);
  chi0reg.SetRealAtm(fRealAtm, fDeex, fEmissionPointCorrection);
  fChi0Regression = &chi0reg;

  // sim chi0 (for testing purposes)
  /*  double chi0, chi0Err, rp, rpErr, t0, t0Err, chi2timing;
  chi0 = eye.GetTelescope(fTelId).GetRecData().GetChiZero();
  chi0Err = 0.01*chi0;
  CombinedChi2(chi0, rp, rpErr, t0, t0Err, PCGFData::eFit, true);
  */

  evt::Event eventCopy1;
  Eye& eyeCopy1 = CopyEye(eye, eventCopy1);

  if (fApLightMethod == eInternal) {
    fApLight.SetStripMode(true);
    const bool haveLight = fApLight.Run(eventCopy1) && eyeCopy1.GetRecData().HasLightFlux();
    if (!haveLight) {
      fProfileFit.UnSetFitterPars();
      return false;
    }
    fEye = &eyeCopy1;
    fPrecise = false;
  }

  double from=fScanStart, to=fScanStop;
  if ((fPrescan?(!Prescan(fScanStep, from, to)):false)
      || !ScanChi0(eyeCopy1, fScanStep, from, to)) {
    fProfileFit.UnSetFitterPars();
    return false;
  }
  const double startchi0 = eyeCopy1.GetTelescope(fTelId).GetRecData().GetChiZero();
  {
    Minou::ParameterDef& chi0Par = GetParameterDefs()[0];
    chi0Par.fValue = startchi0;
    chi0Par.fStep = 0.5*fScanStep;
  }

  // fChi0Regression = &chi0reg; //changed in scanchi0 -- this should be refactored

  double chi0 = startchi0;
  double chi0Err = fScanStep;

  double rp = 0;
  double rpErr = 0;
  double t0 = 0;
  double t0Err = 0;

  evt::Event eventCopy2;
  Eye& eyeCopy2 = CopyEye(eye, eventCopy2);

  if (fApLightMethod == eInternal) {
    CombinedChi2(chi0, rp, rpErr, t0, t0Err, PCGFData::eScan, false);
    FillParams(eyeCopy2, chi0, rp, t0);
    if (!AdjustGeometry(eyeCopy2)) {
      fProfileFit.UnSetFitterPars();
      return false;
    }
    fApLight.SetStripMode(false);
    const bool haveLight = fApLight.Run(eventCopy2) && eyeCopy2.GetRecData().HasLightFlux();
    if (!haveLight) {
      fProfileFit.UnSetFitterPars();
      return false;
    }
    fEye = &eyeCopy2;
    fPrecise = true;
  }

  if (!fScanOnly) {
    Minou::Minimizer<PCGFitter> min(*this);
    if (min.Migrad()) {
      ERROR ("MIGRAD failed.");
      fProfileFit.UnSetFitterPars();
      return false;
    }

    const auto result = min.GetParameterAndError(0);

    chi0 = result.first;
    chi0Err = result.second;
  } else {
    if (!ScanChi0(eyeCopy2, fScanStep/10., startchi0-3*fScanStep, startchi0+3*fScanStep, 1)) {
      fProfileFit.UnSetFitterPars();
      return false;
    }
    chi0 = eyeCopy2.GetTelescope(fTelId).GetRecData().GetChiZero();
    chi0Err = eyeCopy2.GetTelescope(fTelId).GetRecData().GetChiZeroError();
  }

  if (std::isnan(chi0Err) || chi0Err > 180*degree /*||
      chi0 < 5*degree || chi0 > 175*degree*/) {
    cerr << "FIT RESULT:\n\tchi0 = (" << chi0 / degree
         << " +/- " << chi0Err / degree << ") degree\n";
    ERROR ("Fit result not regarded sane.");
    fProfileFit.UnSetFitterPars();
    return false;
  }

  CombinedChi2(chi0, rp, rpErr, t0, t0Err, PCGFData::eFinal, true);

  ostringstream info;
  info << "FIT RESULT:\n\tchi0 = (" << chi0 / degree
       << " +/- " << chi0Err / degree << ") degree\n"
       << "\trp = ( " << rp / km
       << " +/- " << rpErr/km << " ) km\n"
       << "\tt0 = ( " << t0 / microsecond
       << " +/- " << t0Err / microsecond << " ) us";
  INFO(info);

  FillParams(eye, chi0, chi0Err, rp, rpErr, t0, t0Err);
  AdjustGeometry(eye);

  // fill the PCGF data
  {    
    CombinedChi2(chi0, rp, rpErr, t0, t0Err, PCGFData::eFinal, true);
    
    EyeRecData& eyeRecWrite = eye.GetRecData();
    eyeRecWrite.SetPCGF(fPCGFData);
  }

  fProfileFit.UnSetFitterPars();

  return true;
}


bool
PCGFitter::Prescan(const double step, double &from, double &to){
  //calling CombinedChi2, so fChi0Regression must be set
  const double bigstep=4*step;
  const double inf = numeric_limits<double>::infinity();

  double rp=0, t0=0, rpErr=0, t0Err=0;
  for (;; from+=bigstep) {
    if (from >= to) {
      PCGF_DEBUGLOG("Pescan failed leave: All profile fits failed.");
      return false;
    }

    if (CombinedChi2(from, rp, rpErr, t0, t0Err, PCGFData::eScan, true) < inf) {
      break;
    }
  }
  //to be sure that we do not miss anything
  from -= bigstep;

  for (;; to-=bigstep) {
    if (from > to) {
      PCGF_DEBUGLOG("Pescan failed leave: All profile fits failed backwards. This should not happen.");
      return false;
    }

    if (CombinedChi2(to, rp, rpErr, t0, t0Err, PCGFData::eScan, true) < inf) {
      break;
    }
  }
  to += bigstep;
  PCGF_DEBUGLOG( "Prescan: chi^2 finite for chi0 from : "<< from << "\t to: " << to);
  if (from >= to) {
    return false;
  }

  return true;
}
  
bool
PCGFitter::ScanChi0(Eye &eye, const double step, double from, double to, int level)
{
  const double inf = numeric_limits<double>::infinity();

  PCGF_DEBUGLOG("ScanChi0 enter");
  //const EyeRecData& eyeRec = eye.GetRecData();
  //ChiZeroRegression chi0reg(eyeRec);
  //chi0reg.SetUseLightFlux(fUseLightFlux);
  //fChi0Regression = &chi0reg;
  //ChiZeroRegression &chi0reg = *fChi0Regression;

  vector<double> chi0s;
  vector<double> chi2s;
  vector<double> rps;
  vector<double> t0s;

  for (double chi0 = from; chi0 <= to; chi0 += step) {
    double rp = 0, rpErr = 0;
    double t0 = 0, t0Err = 0;
    double chi2 = CombinedChi2(chi0, rp, rpErr, t0, t0Err, PCGFData::eScan, true);

    chi0s.push_back(chi0);
    chi2s.push_back(chi2);
    rps.push_back(rp);
    t0s.push_back(t0);
  }

  unsigned int min_i = TMath::LocMin(chi2s.size(), &chi2s[0]);
  
  PCGF_DEBUGLOG( "ScanChi0: min i:" << min_i << "\t size: " << chi2s.size());
  PCGF_DEBUGLOG( "ScanChi0: min chi2: " << chi2s[min_i]);

  if (!(chi2s[min_i] < inf)) {
    PCGF_DEBUGLOG("ScanChi0 failed leave: All profile fits failed.");
    return false;
  }

  if (level > 0 && level < 10) { //we are in the second scan phase (do at most 10 re-scans)
    //if we are at the scan range border, do another scan 
    if (min_i == 0) {
      return ScanChi0(eye, fScanStep/10., chi0s[min_i]-3*fScanStep, chi0s[min_i]+4*fScanStep/10., level+1);
    } else if (min_i == chi2s.size()-1) {
      return ScanChi0(eye, fScanStep/10., chi0s[min_i]-4*fScanStep/10., chi0s[min_i]+3*fScanStep, level+1);
    }
  }
/*  else if (level == 4 && (min_i == 0 || min_i == chi2s.size()-1)) {
    PCGF_DEBUGLOG("ScanChi0 failed leave: No smooth minimum after 3 re-scans.");
    return false;
  }
*/

  /*const double tooMuch = numeric_limits<double>::max(); //but not infinity! //TODO: ask Vladimir why
  //VN : imagine example [..66,55,33,inf..] -> [..66,55,tooMuch,inf..] -> 55 is OK
  //     if tooMuch=inf :[..66,55,33,inf..] -> [..66,55,inf,inf..] -> [..66,inf,inf,inf..] etc.
  while (!((min_i > 0) ? chi2s[min_i-1] : inf < inf &&
         (min_i < chi2s.size()-1) ? chi2s[min_i+1] : inf < inf) &&
         (chi2s[min_i] != tooMuch)) { //What about the margins? - VN : true, fixed
    chi2s[min_i] = tooMuch;
    min_i = TMath::LocMin(chi2s.size(), &chi2s[0]);
  }
  if (chi2s[min_i] == tooMuch) {
    PCGF_DEBUGLOG("ScanChi0 failed leave: No smooth minimum.");
    return false;
  }*/

  //VN: different estimate of smooth minimum
  const double tooMuch = numeric_limits<double>::max();
  vector<double> chi2s_copy(chi2s);
  chi2s_copy[0] = tooMuch;
  chi2s_copy[chi2s_copy.size()-1] = tooMuch;
  min_i = TMath::LocMin(chi2s_copy.size(), &chi2s_copy[0]);

  while (level>0 && chi2s_copy[min_i] < tooMuch) {
    vector<double> chi2s_trial(chi2s_copy);
    chi2s_trial[min_i] = tooMuch;
    int min_i_trial1 = TMath::LocMin(chi2s_trial.size(), &chi2s_trial[0]);
    chi2s_trial[min_i_trial1] = tooMuch;
    int min_i_trial2 = TMath::LocMin(chi2s_trial.size(), &chi2s_trial[0]);
    chi2s_trial[min_i_trial2] = tooMuch;
    if ((chi2s_trial[min_i-1] == tooMuch) && (chi2s_trial[min_i+1] == tooMuch)) {
      break;
    } else {
      chi2s_copy[min_i] = tooMuch;
      min_i = TMath::LocMin(chi2s_copy.size(), &chi2s_copy[0]);
    }
  }

  if (!(chi2s_copy[min_i] < tooMuch)) {
    PCGF_DEBUGLOG("ScanChi0 failed leave: No smooth minimum.");
    return false;
  }

  //estimate chi0 error
  double xmin = chi0s[min_i];
  double ymin = chi2s[min_i];
  vector<double> x;
  vector<double> y;
  vector<double> yerr;
  x.push_back(0.);
  y.push_back(0.);
  yerr.push_back(1.);
  for (unsigned int i = 1; i < 5; ++i) {
    int j = min_i-i;
    unsigned int k = min_i+i;
    if ((j > 0) && std::isfinite(chi2s[j]) && (chi2s[j] != tooMuch)) {
      x.push_back(abs(chi0s[j]-xmin));
      y.push_back(((chi2s[j]-ymin)>0.) ? sqrt(chi2s[j]-ymin) : 0.);
      yerr.push_back(1.);
    //cout<<x[x.size()-1] << y[y.size()-1]<<endl;
    }
    if ((k < chi0s.size()) && std::isfinite(chi2s[k]) && (chi2s[k] != tooMuch)) {
      x.push_back(abs(chi0s[k]-xmin));
      y.push_back(((chi2s[k]-ymin)>0.) ? sqrt(chi2s[k]-ymin) : 0.);
      yerr.push_back(1.);
    //cout<<x[x.size()-1] << y[y.size()-1]<<endl;
    }
  }
  double a = 0.;
  double b = 0.;
  double chi2 = 0.;
  otoa::LinearFit(x, y, yerr, a, b, chi2);
  double chi0err = (b!=0) ? (1./b) : 90*deg;
  //end of chi0 error estimation

  FillParams(eye, chi0s[min_i], chi0err, rps[min_i], 0.1*kilometer, t0s[min_i], 0.1*microsecond);

  cout << "calculating result geometry\n";
  AdjustGeometry(eye);
  {
    const evt::ShowerFRecData& recShower = eye.GetRecData().GetFRecShower();
    const Point& corePos = recShower.GetCorePosition();
    const CoordinateSystemPtr coreCS = fwk::LocalCoordinateSystem::Create(corePos);
    const Vector& axis = recShower.GetAxis();
    cout << "chi0: " << chi0s[min_i]/degree << " deg  "
            "rp: " << rps[min_i] / meter << " m  "
            " -> "
            "theta: " << axis.GetTheta(coreCS)/degree << " deg  "
            "phi: " << axis.GetPhi(coreCS)/degree << " deg"
         << endl;
  }

  PCGF_DEBUGLOG("ScanChi0 sucessful leave");
  return true;
}


Eye&
PCGFitter::CopyEye(const Eye& eye, evt::Event& eventCopy)
  const
{
  eventCopy.MakeFEvent();
  eventCopy.GetFEvent().MakeEye(eye.GetId());
  eventCopy.GetFEvent().GetEye(eye.GetId()) = eye;
  return eventCopy.GetFEvent().GetEye(eye.GetId());
}


void
PCGFitter::FillParams(Eye& eyeCopy,
                      const double chi0, const double chi0Err,
                      const double rp, const double rpErr,
                      const double t0, const double t0Err)
  const
{
  if (eyeCopy.HasTelescope(fTelId)) {
    fevt::Telescope& tel = eyeCopy.GetTelescope(fTelId);
    if (!tel.HasRecData())
        tel.MakeRecData();
    fevt::TelescopeRecData& telRecData = tel.GetRecData();
    telRecData.SetChiZero(chi0,chi0Err);
    telRecData.SetTZero(t0,t0Err);
    telRecData.SetRp(rp,rpErr);
    telRecData.SetTimeFitChiSquare(0,1);
    telRecData.SetTimeFitCorrelations(0,0,0);
  }
}


void
PCGFitter::FillParams(Eye& eyeCopy,
                      const double chi0,
                      const double rp,
                      const double t0)
  const
{
  FillParams(eyeCopy, chi0, 0.5*degree, rp, 0.1*km, t0, 0.1*microsecond);
}


bool
PCGFitter::AdjustGeometry(Eye& eyeCopy)
  const
{
  
  EyeRecData& eyeRecData = eyeCopy.GetRecData();
  const fdet::FDetector& detFD = det::Detector::GetInstance().GetFDetector();
  const fdet::Eye& detEye = detFD.GetEye(eyeCopy);
  const Point eyePos = detEye.GetPosition();

  const utl::TimeStamp eyeStamp = eyeCopy.GetHeader().GetTimeStamp();

  if (!eyeCopy.HasTelescope(fTelId) || !eyeCopy.GetTelescope(fTelId).HasRecData()) {
    return false;
  }
  TelescopeRecData& telRecData = eyeCopy.GetTelescope(fTelId).GetRecData();

  double t0 = telRecData.GetTZero();
  const double chi0 = telRecData.GetChiZero();
  const double rp = telRecData.GetRp();
  const double t0Err = telRecData.GetTZeroError();
  const double chi0Err = telRecData.GetChiZeroError();
  const double rpErr = telRecData.GetRpError();

  const Point telPos = detEye.GetTelescope(fTelId).GetPosition();
  const CoordinateSystemPtr telCS = fwk::LocalCoordinateSystem::Create(telPos);
  const double sdpTheta = telRecData.GetSDP().GetTheta(telCS);
  const double sdpPhi = telRecData.GetSDP().GetPhi(telCS);

  const Vector sdp = Vector(1, sdpTheta, sdpPhi, telCS, Vector::kSpherical);
  const Vector vertical(0,0,1, telCS);
  const Vector horizontalInSDP = Normalized(Cross(sdp, vertical));

  const Transformation rot = Transformation::Rotation(-chi0, sdp, telCS);
  Vector axis = Normalized(rot*horizontalInSDP);

  const double coreDistance = rp / sin(chi0);
  const Vector coreTelVec = coreDistance * horizontalInSDP;
  const Point core = telPos + coreTelVec;
  const Point P_in_tel = telPos + Normalized(Cross(sdp, axis))*rp;

  //VN: set telescope data to deal with HECO; also set eye data and recalculate chi_i
  bool fillEye = true;
  for (fdet::Eye::TelescopeIterator detTelIt = detEye.TelescopesBegin();
       detTelIt != detEye.TelescopesEnd(); ++detTelIt)
  {

    const int curtelId = detTelIt->GetId();

    const Point curtelPos = detTelIt->GetPosition();
    const CoordinateSystemPtr curtelCS = fwk::LocalCoordinateSystem::Create(curtelPos);

    const Point core_in_curtel = core-axis*core.GetZ(curtelCS)/cos(axis.GetTheta(curtelCS));

    const Vector curtel_to_core = core_in_curtel-curtelPos;
    double chi0_in_curtel = Angle(curtel_to_core, axis);
    if (axis.GetTheta(curtelCS)>kPi/2.) {
      chi0_in_curtel = -(kPi-chi0_in_curtel);
    }
    double rp_in_curtel = curtel_to_core.GetMag()*sin(chi0_in_curtel);
    if (axis.GetTheta(curtelCS)>kPi/2.) {
      rp_in_curtel *= -1.;
    }
    const utl::AxialVector SDP_in_curtel = Normalized(Cross(axis, curtel_to_core));
    const Point P_in_curtel = curtelPos + Normalized(Cross(SDP_in_curtel, axis))*rp_in_curtel;
    const Vector Pt_to_Pct = P_in_curtel - P_in_tel;
    const double t0_in_curtel = t0 + ((Pt_to_Pct*axis > 0.) ? -1. : 1.)*Pt_to_Pct.GetMag()/kSpeedOfLight;

    if (eyeCopy.HasTelescope(curtelId)) {
      fevt::Telescope& tel = eyeCopy.GetTelescope(curtelId);

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

      fevt::TelescopeRecData& telRecData = tel.GetRecData();
      telRecData.SetSDP(SDP_in_curtel);
      telRecData.SetChiZero(chi0_in_curtel, chi0Err);
      telRecData.SetTZero(t0_in_curtel, t0Err);
      telRecData.SetRp(rp_in_curtel, rpErr);
      telRecData.SetTimeFitChiSquare(0,1);
      telRecData.SetTimeFitCorrelations(0,0,0);
    }

    if (fillEye) {
      //choose tel with the same position as eye (usually, more of them are available)
      if ((curtelPos - eyePos).GetMag() < 1e-3) {
        eyeRecData.SetChiZero(chi0_in_curtel, chi0Err);
        eyeRecData.SetTZero(t0_in_curtel, t0Err);
        eyeRecData.SetRp(rp_in_curtel, rpErr);
        eyeRecData.SetSDP(SDP_in_curtel);

        evt::ShowerFRecData& frecShower = eyeRecData.GetFRecShower();
        frecShower.SetAxis(axis);
        frecShower.SetCorePosition(core_in_curtel);
        try {
          double tCore = t0_in_curtel;
          tCore -= (rp_in_curtel / tan(chi0_in_curtel)) / kSpeedOfLight;
          utl::TimeStamp timeAtCore = eyeStamp - TimeInterval(1000*100*ns) + TimeInterval(tCore);
          frecShower.SetCoreTime(timeAtCore, 0); //, rp_in_tel, chi0_in_tel, t0_in_tel);
        } catch (OutOfBoundException& e) {
          ERROR(e.GetMessage());
          return false;
        }
        fillEye = false;
      }
    }
  }

  for (fevt::EyeRecData::PixelIterator pixelIter = eyeRecData.PulsedPixelsBegin();
      pixelIter != eyeRecData.PulsedPixelsEnd(); ++pixelIter) {

    unsigned int telId = pixelIter->GetTelescopeId();

    fevt::TelescopeRecData& telRecData = eyeCopy.GetTelescope(telId).GetRecData();
    const Point telPos = det::Detector::GetInstance().GetFDetector().GetEye(eyeCopy).GetTelescope(telId).GetPosition();
    const CoordinateSystemPtr telCS = fwk::LocalCoordinateSystem::Create(telPos);
    const double sdpTheta = telRecData.GetSDP().GetTheta(telCS);
    const double sdpPhi = telRecData.GetSDP().GetPhi(telCS);

    const Vector sdp = Vector(1, sdpTheta, sdpPhi, telCS, Vector::kSpherical);
    const Vector vertical(0,0,1, telCS);
    const Vector horizontalInSDP = Normalized(Cross(sdp, vertical));

    const Vector pixelDirection = detFD.GetPixel(*pixelIter).GetDirection();
    const Vector pixelDirInSDP = pixelDirection - (sdp * pixelDirection)*sdp;
    const double pixelChi = Angle(pixelDirInSDP, horizontalInSDP);

    pixelIter->GetRecData().SetChi_i(pixelChi);
  }

  return true;
}

}