G/ChiZeroRegression.cc

Go to the documentation of this file.

#include "ChiZeroRegression.h"

#include <utl/Point.h>
#include <utl/Vector.h>
#include <utl/CoordinateSystemPtr.h>
#include <utl/Transformation.h>
#include <utl/UTMPoint.h>
#include <fwk/LocalCoordinateSystem.h>

#include <utl/PhysicalConstants.h>
#include <utl/PhysicalFunctions.h>
#include <fevt/Pixel.h>
#include <fevt/PixelRecData.h>
#include <fevt/Telescope.h>
#include <fevt/TelescopeRecData.h>
//#include <fevt/EyeRecData.h>

#include <det/Detector.h>
#include <fdet/FDetector.h>
#include <fdet/Eye.h>
#include <fdet/Telescope.h>
#include <fdet/Pixel.h>
#include <utl/CoordinateSystemPtr.h>
#include <fdet/FTimeFitModel.h>
#include <atm/Atmosphere.h>
#include <atm/InclinedAtmosphericProfile.h>
#include <atm/ProfileResult.h>
#include <utl/ReferenceEllipsoid.h>
#include <TMinuit.h>
#include <TMath.h>

#include <limits>

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

namespace FdProfileConstrainedGeometryFitPG {

int ChiZeroRegression::fEyeId = 0;
int ChiZeroRegression::fTelId = 0;
double ChiZeroRegression::fXmax = 0;
double ChiZeroRegression::fsdpTheta = 0;
double ChiZeroRegression::fsdpPhi = 0;
std::vector<double> ChiZeroRegression::times;
std::vector<double> ChiZeroRegression::timeErrs;
std::vector<double> ChiZeroRegression::chii;
std::vector<int> ChiZeroRegression::tels;
std::vector<int> ChiZeroRegression::pixids;
bool ChiZeroRegression::fEmissionPointCorrection = false;

void
ChiZeroRegression::SetRealAtm(bool realAtm, bool deex, bool emissionPointCorrection) {
  fRealAtm = realAtm;
  fDeex = deex;
  fdet::FTimeFitModel::GetInstance().SetModel((fRealAtm && fDeex) ? (fdet::FTimeFitModel::eRealisticAtm) : 
        ((fRealAtm && !fDeex) ? (fdet::FTimeFitModel::eRealisticAtmNoDeex) : (fdet::FTimeFitModel::eVacuumAtm)));
  fEmissionPointCorrection = emissionPointCorrection;
}

unsigned int
ChiZeroRegression::PrepParams(const fevt::Eye& eye)
{
  using namespace fevt;

  const EyeRecData& eyeRec = eye.GetRecData();
  const fdet::FDetector& detFD = det::Detector::GetInstance().GetFDetector();
  const unsigned int lastTel = detFD.GetEye(eyeRec.GetEyeId()).GetLastTelescopeId();
  std::vector<unsigned int> nPix(lastTel, 0);

  for (EyeRecData::ConstPixelIterator pix = eyeRec.TimeFitPixelsBegin(),
         end = eyeRec.TimeFitPixelsEnd(); pix != end; ++pix) {
    nPix[pix->GetTelescopeId()-1]++;
  }
  unsigned int maxPix = 0;
  for (unsigned int i=0; i<nPix.size(); i++) {
    if (nPix[i] > maxPix) {
      fTelId = i+1;
      maxPix = nPix[i];
    }
  }

  fEyeId = eyeRec.GetEyeId();

  chii.resize(0);
  times.resize(0);
  timeErrs.resize(0);
  tels.resize(0);
  pixids.resize(0);

  for (EyeRecData::ConstPixelIterator pix = eyeRec.TimeFitPixelsBegin(),
         end = eyeRec.TimeFitPixelsEnd(); pix != end; ++pix) {
    const PixelRecData& pixRec = pix->GetRecData();
    chii.push_back(pixRec.GetChi_i());
    times.push_back(pixRec.GetT_i().GetNanoSecond());
    timeErrs.push_back(pixRec.GetT_iError().GetInterval());
    tels.push_back(pix->GetTelescopeId());
    pixids.push_back(pix->GetId());
  }

  const utl::Point telPos = detFD.GetEye(fEyeId).GetTelescope(fTelId).GetPosition();
  utl::CoordinateSystemPtr telCS = fwk::LocalCoordinateSystem::Create(telPos);
  if (eye.HasTelescope(fTelId) && eye.GetTelescope(fTelId).HasRecData()) {
    fsdpTheta = eye.GetTelescope(fTelId).GetRecData().GetSDP().GetTheta(telCS);
    fsdpPhi = eye.GetTelescope(fTelId).GetRecData().GetSDP().GetPhi(telCS);
  }

  return fTelId;
}

void
ChiZeroRegression::operator()(const double chi0,
                              double& chi2,
                              double& rp,
                              double& t0,
                              double Xmax)
  const
{
  std::vector<double> xs;//(chii.size());
  std::vector<double> ys;
  std::vector<double> yerrs;

  double rpErr = 0.;
  double t0Err = 0.;
  for (unsigned int i = 0, n = chii.size(); i < n; ++i) {
    if (tels[i] == fTelId) {
      xs.push_back(tan(0.5*(chi0 - chii[i])));
      ys.push_back(times[i]);
      yerrs.push_back(timeErrs[i]);
    }
  }

  if (!fUseLightFlux) {
    otoa::LinearFit(xs, ys, yerrs, t0, rp, chi2);
  } else {
    LinearFitErrors(xs, ys, yerrs, t0, t0Err, rp, rpErr, chi2);
  }
    
  rp *= utl::kSpeedOfLight;
  rpErr *= utl::kSpeedOfLight;

  fXmax = Xmax;
//  if (fRealAtm)
    MinuitFitErrors(chi0, t0, t0Err, rp, rpErr, chi2);
}

void
ChiZeroRegression::operator()(const double chi0,
                              double& chi2,
                              double& rp, double &rpErr,
                              double& t0, double &t0Err,
                              double Xmax)
  const
{
  std::vector<double> xs;//(chii.size());
  std::vector<double> ys;
  std::vector<double> yerrs;

  for (unsigned int i = 0, n = chii.size(); i < n; ++i) {
    if (tels[i] == fTelId) {
      xs.push_back(tan(0.5*(chi0 - chii[i])));
      ys.push_back(times[i]);
      yerrs.push_back(timeErrs[i]);
    }
  }

  LinearFitErrors(xs, ys, yerrs, t0, t0Err, rp, rpErr, chi2);
  rp *= utl::kSpeedOfLight;
  rpErr *= utl::kSpeedOfLight;

  fXmax = Xmax;
//  if (fRealAtm)
    MinuitFitErrors(chi0, t0, t0Err, rp, rpErr, chi2);
}

void 
ChiZeroRegression::MinuitFitFunc(int& /*npar*/, double* /*gin*/, double& f, double* par, int /*iflag*/)
{
  using namespace utl;

  const double dirErr = 0.35*deg;

  double& Chi0 = par[0];
  double& T0 = par[1];
  double& Rp = par[2];
  fdet::FTimeFitModel& timeFitModel = fdet::FTimeFitModel::GetInstance();

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

  const Point telPos = detFD.GetEye(fEyeId).GetTelescope(fTelId).GetPosition();
  utl::CoordinateSystemPtr telCS = fwk::LocalCoordinateSystem::Create(telPos);
  const Vector sdp = Vector(1, fsdpTheta, fsdpPhi, 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;

  //correction for shower LDF
  try {
    const atm::Atmosphere& atmo = det::Detector::GetInstance().GetAtmosphere();
    atmo.InitSlantProfileModel(core, axis, 10*g/cm2);
    const atm::ProfileResult& tempProfile = atmo.EvaluateTemperatureVsHeight();
    const atm::ProfileResult& pressureProfile = atmo.EvaluatePressureVsHeight();
    const atm::ProfileResult& slantDepthProfile = atmo.EvaluateSlantDepthVsDistance();
    const ReferenceEllipsoid& wgs84 = ReferenceEllipsoid::GetWGS84();

  double logLikeSum = 0.;
  double logLikeSumSDP = 0.;
  for (int i = 0, n = chii.size(); i < n; ++i) {

    const fdet::Telescope& curTel = detFD.GetEye(fEyeId).GetTelescope(tels[i]);
    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);
    }
    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;

    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(pixids[i]).GetDirection();
    const Vector pixelDirInSDP = pixelDirection - (SDP_in_curtel * pixelDirection)*SDP_in_curtel;
    const double pixelChi = Angle(pixelDirInSDP, horizontalInSDP_in_curtel);

    const double angleWithSDP = abs(kPiOnTwo - Angle(SDP_in_curtel, pixelDirection));

    const Point PointOnShower = P_in_curtel+axis*rp_in_curtel/tan(chi0_in_curtel-pixelChi);
    const double cosLocalZenith =
          axis.GetCosTheta(fwk::LocalCoordinateSystem::Create(UTMPoint(PointOnShower, wgs84)));
    const double height = wgs84.PointToLatitudeLongitudeHeight(PointOnShower).get<2>();
    const double pressure = pressureProfile.Y(height);
    const double temperature = tempProfile.Y(height);
    const double rMoliere = MoliereRadius(temperature, pressure, cosLocalZenith);
    const Vector core_to_PoS = PointOnShower-core;
    const double age = ShowerAge(slantDepthProfile.Y( ((core_to_PoS*axis > 0.) ? -1. : 1.) * core_to_PoS.GetMag() ),fXmax);
    // parameters based on 3D CHERENKOV CORSIKA simulations at low energies (10^15-10^17 eV)
    const double rcorr = InverseGoraCDF(0.9, age)*0.2954;
    const double rMcor = rMoliere*rcorr;
    double t_corr = (-rMcor/kSpeedOfLight*tan((chi0_in_curtel-pixelChi)/2.))+5.923*m/kSpeedOfLight;

    double ti_exp = timeFitModel.GetTimeAtAperture(t0_in_curtel, rp_in_curtel, chi0_in_curtel, pixelChi, 
                                   SDP_in_curtel.GetTheta(curtelCS), fEyeId, tels[i]) + (fEmissionPointCorrection ? t_corr : 0.);
    logLikeSum += (-0.5*(2*kPi) - log(timeErrs[i]) - 0.5*pow( (times[i] - ti_exp)/timeErrs[i] , 2));
    //SDP part
    logLikeSumSDP += (-0.5*(2*kPi) - log(dirErr) - 0.5*pow(angleWithSDP/dirErr, 2));
  }
  logLikeSum += logLikeSumSDP;//*chii.size();
  if (logLikeSum != 0) {
    f = -2.*logLikeSum;
  } else {
    f = std::numeric_limits<double>::infinity();
  }

  } catch (const std::exception&) {
    f = std::numeric_limits<double>::infinity();
    return;
  }
}

void
ChiZeroRegression::MinuitFitErrors(const double Chi0,
                                   double& T0, double& T0err,
                                   double& Rp, double& Rperr,
                                   double& chi2) const
{
  const int npar = 3;
  TMinuit theMinuit(npar);
  theMinuit.SetPrintLevel(-1);
  theMinuit.SetFCN(ChiZeroRegression::MinuitFitFunc);
  double arglist[npar];
  int ierflag=0;
  arglist[0] = 1;
  arglist[1] = 1;
  theMinuit.mnexcm("SET ERR", arglist,1,ierflag);
  if (ierflag) ERROR("Minuit SET ERR failed");
  static double vstart[npar];
  static double step[npar];
  vstart[0] = Chi0;
  vstart[1] = T0;
  vstart[2] = Rp;
  step[0] = 0;
  step[1] = 10;
  step[2] = 10;
  theMinuit.mnparm(0, "Chi0", vstart[0], step[0], 0,0,ierflag);
  theMinuit.FixParameter(0);
  theMinuit.mnparm(1, "T0", vstart[1], step[1], vstart[1]-300, vstart[1]+300.,ierflag);  //hardcoded boundaries, should be optional in the future?
  theMinuit.mnparm(2, "Rp", vstart[2], step[2], vstart[2]-300., vstart[2]+300.,ierflag);
  arglist[0] = 3000; //maxcalls
  arglist[1] = 1;    //tolerance
  theMinuit.mnexcm("MIGRAD", arglist, 2, ierflag);
  if (ierflag) {
    ERROR("Minuit MIGRAD failed");
  }  
  theMinuit.GetParameter(1, T0, T0err);
  theMinuit.GetParameter(2, Rp, Rperr);
  chi2 = theMinuit.fAmin;
}

void
ChiZeroRegression::LinearFitErrors(const std::vector<double>& x,
                                   const std::vector<double>& y,
                                   const std::vector<double>& ey,
                                   double& a0, double& a0err,
                                   double& a1, double& a1err,
                                   double& chi2) const
{
  double sw = 0;
  double sx = 0;
  double sxx = 0;
  double sy = 0;
  double sxy = 0;
  double syy = 0;

  for (unsigned int i = 0; i<x.size(); ++i) {
    if (ey[i] > 0) {
      const double wi = 1.0 / (ey[i]*ey[i]);
      sw  += wi;
      sx  += wi * x[i];
      sxx += wi * x[i] * x[i];
      sy  += wi * y[i];
      sxy += wi * x[i] * y[i];
      syy += wi * y[i] * y[i];
    }
  }

  const double dinv = 1.0/(sw*sxx - sx*sx);

  a0 = (sxx*sy - sx * sxy) * dinv;
  a1 = (-sx*sy + sw * sxy) * dinv;
  a0err = sqrt (sxx * dinv);
  a1err = sqrt (sw * dinv);

  if (!fUseLightFlux) {
    chi2 = syy - a0*sy - a1*sxy;
  } else {
    double logLikeSum = 0; //differs only by constant from chi2, but put in
    for (unsigned int i = 0; i<x.size(); ++i) {
      logLikeSum += (-0.5*log(ey[i]*ey[i]) - 0.5*pow( (y[i]-(a0+a1*x[i]))/ey[i] , 2));
    }
    chi2 = -2.*logLikeSum;
  }
}

} //namespace