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 <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 FdProfileConstrainedGeometryFit {

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

  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("delZeroLightFlux").GetData(fDelZeroFlux);
  topBranch.GetChild("TimeFitModel").GetData(fRealAtm);
  topBranch.GetChild("TimeFitDeexcitation").GetData(fDeex);
  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 TimeFitModel is " << (fRealAtm?"RealisticAtm":"VacuumAtm") << "\n"
          "\t TimeFitDeexcitation is " << (fDeex?"on":"off") << "\n";
  INFO(info);

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

  fProfileFit.Init();

  //  fChi0Regression->
  fProfileFit.SetUseLightFlux(fUseLightFlux);
}


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

  const EyeRecData& eyeRec = eye.GetRecData();

  for (EyeRecData::ConstPixelIterator pix = eyeRec.TimeFitPixelsBegin(),
       end = eyeRec.TimeFitPixelsEnd(); pix != end; ++pix) {
    const PixelRecData& pixRec = pix->GetRecData();
    const double chii = pixRec.GetChi_i();
    //cout << "chii " << chii << "pixel " << (*pix).GetId() << "\n";
    if (sin(chi0 - 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;
  const double chi2 = CombinedChi2(chi0, rp, t0, PCGFData::eFit, true);
  return chi2;
}


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

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

  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, 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]);
        }
      }
    }
  }

  const double chi2profile = fProfileFit(eyeCopy);
  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';

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

  // {
  //   const double startchi0=ScanChi0(eye, fScanStep,fScanStart,fScanStop)?
  //     eye.GetRecData().GetChiZero() : 90*degree;
  //   SetInitialValueStep(0, startchi0, 3.0*degree);
  // }

  const EyeRecData& eyeRec = eye.GetRecData();
  ChiZeroRegression chi0reg(eyeRec);
  chi0reg.SetUseLightFlux(fUseLightFlux);
  chi0reg.SetRealAtm(fRealAtm, fDeex);
  fChi0Regression = &chi0reg;

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

  double from = fScanStart;
  double to = fScanStop;
  if ((fPrescan && !Prescan(fScanStep, from, to)) || !ScanChi0(eye, fScanStep, from, to)) {
    fProfileFit.UnSetFitterPars();
    return false;
  }
  const double startchi0 = eye.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 chi2timing = 0;
  double rp = 0;
  double rpErr = 0;
  double t0 = 0;
  double t0Err = 0;

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

  if (fApLightMethod == eInternal) {
    chi0reg(chi0, chi2timing, rp, t0);
    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;
  }

  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(eye, 0.1*fScanStep, startchi0 - 3*fScanStep, startchi0 + 3*fScanStep)) {
      fProfileFit.UnSetFitterPars();
      return false;
    }
    chi0 = eye.GetRecData().GetChiZero();
    chi0Err = eye.GetRecData().GetChiZeroError();
  }

  if (std::isnan(chi0Err) || chi0Err > 90*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;
  }

  chi0reg(chi0, chi2timing, rp, rpErr, t0, t0Err);

  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, t0, 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;
  double t0 = 0;
  for ( ; ; from += bigstep) {
    if (from >= to) {
      PCGF_DEBUGLOG("Pescan failed leave: All profile fits failed.");
      return false;
    }

    if (CombinedChi2(from, rp, t0, 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, t0, 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)
{
  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;
    double t0 = 0;
    double chi2 = CombinedChi2(chi0, rp, t0, PCGFData::eScan, true);

    chi0s.push_back(chi0);
    chi2s.push_back(chi2);// > 0.1 ? chi2 : inf); // WARNING: what is this check good for? - removed
    rps.push_back(rp);
    t0s.push_back(t0);
  }

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

  /*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();
  while (chi2s[min_i] < tooMuch) {
    vector<double> chi2s_trial(chi2s);
    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[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;
  }
  /* check perhaps not needed any more: - VN : not needed + eye not filled yet => drop off showers
  if (!(chi2s[min_i] < numeric_limits<double>::infinity()) ||
      eye.GetRecData().GetChiZero() <= (from+2*step) ||
      eye.GetRecData().GetChiZero() >= (to-2*step)) {
    PCGF_DEBUGLOG("ScanChi0 failed leave: No minimum in range-margin");
    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(std::abs(chi0s[j] - xmin));
      y.push_back(((chi2s[j] - ymin) > 0) ? sqrt(chi2s[j] - ymin) : 0.);
      yerr.push_back(1.);
    }
    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.);
    }
  }
  double a = 0.;
  double b = 0.;
  double chi2 = 0.;
  otoa::LinearFit(x, y, yerr, a, b, chi2);
  double chi0err = (b != 0) ? (1./b) : 90*degree;
  //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
{
  EyeRecData& eyeReco = eyeCopy.GetRecData();
  eyeReco.SetChiZero(chi0, chi0Err);
  eyeReco.SetTZero(t0, t0Err);
  eyeReco.SetRp(rp, rpErr );

  //VN: set telescope data to show real geometry in ti-chii plot in EventBrowser
  for (fevt::Eye::TelescopeIterator telIt = eyeCopy.TelescopesBegin(ComponentSelector::eHasData);
       telIt != eyeCopy.TelescopesEnd(ComponentSelector::eHasData); ++telIt) {
    if (!telIt->HasRecData())
      telIt->MakeRecData();
    fevt::TelescopeRecData& telRecData = telIt->GetRecData();
    telRecData.SetChiZero(chi0,chi0Err);
    telRecData.SetTZero(t0,t0Err);
    telRecData.SetRp(rp,rpErr);
    telRecData.SetTimeFitChiSquare(0,0);
    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 double t0 = eyeRecData.GetTZero();
  const double chi0 = eyeRecData.GetChiZero();
  const double rp = eyeRecData.GetRp();

  const fdet::Eye& detEye =
    det::Detector::GetInstance().GetFDetector().GetEye(eyeCopy);
  const CoordinateSystemPtr eyeCS = detEye.GetEyeCoordinateSystem();
  const double sdpTheta = eyeRecData.GetSDP().GetTheta(eyeCS);
  const double sdpPhi = eyeRecData.GetSDP().GetPhi(eyeCS);

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

  const Transformation rot = Transformation::Rotation(-chi0, sdp, eyeCS);
  const Vector axis = Normalized(rot*horizontalInSDP);
  const double coreDistance = rp / sin(chi0);
  const Vector coreEyeVec = coreDistance * horizontalInSDP;
  const Point& eyepos = detEye.GetPosition();
  const Point core = eyepos + coreEyeVec;
  evt::ShowerFRecData& recShower = eyeRecData.GetFRecShower();
  recShower.SetAxis(axis);
  recShower.SetCorePosition(core);
  try {
    recShower.SetCoreTime(eyeCopy.GetHeader().GetTimeStamp() -
                          TimeInterval(1000*100*ns), rp, chi0, t0);
  } catch (OutOfBoundException& e) {
    ERROR(e.GetMessage());
    return false;
  }
  return true;
}

}