ChiZeroRegression.cc

Go to the documentation of this file.

#include "ChiZeroRegression.h"

#include <utl/PhysicalConstants.h>
#include <fevt/Pixel.h>
#include <fevt/PixelRecData.h>

#include <det/Detector.h>
#include <fdet/FDetector.h>
#include <fdet/Eye.h>
#include <utl/CoordinateSystemPtr.h>
#include <fdet/FTimeFitModel.h>
#include <TMinuit.h>

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

namespace FdProfileConstrainedGeometryFit {

int ChiZeroRegression::fEyeId = 0;
int ChiZeroRegression::fTelId = 0;
double ChiZeroRegression::fsdpTheta = 0;
std::vector<double> ChiZeroRegression::times;
std::vector<double> ChiZeroRegression::timeErrs;
std::vector<double> ChiZeroRegression::chii;

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

void
ChiZeroRegression::PrepParams(const fevt::EyeRecData& eyeRec)
{
  using namespace fevt;

  fEyeId = eyeRec.GetEyeId();
  fTelId = eyeRec.TimeFitPixelsBegin()->GetTelescopeId();
  chii.resize(0);
  times.resize(0);
  timeErrs.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());
    timeErrs.push_back(pixRec.GetT_iError());
  }

  utl::CoordinateSystemPtr eyeCS = det::Detector::GetInstance().GetFDetector().GetEye(fEyeId).GetEyeCoordinateSystem();
  fsdpTheta = eyeRec.GetSDP().GetTheta(eyeCS);
}

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

  double rpErr = 0.;
  double t0Err = 0.;
  for (unsigned int i = 0, n = chii.size(); i < n; ++i)
    xs[i] = tan(0.5*(chi0 - chii[i]));
  if (!fUseLightFlux) {
    otoa::LinearFit(xs, times, timeErrs, t0, rp, chi2);
  } else {
    LinearFitErrors(xs, times, timeErrs, t0, t0Err, rp, rpErr, chi2);
  }
    
  rp *= utl::kSpeedOfLight;
  rpErr *= utl::kSpeedOfLight;

  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)
  const
{
  std::vector<double> xs(chii.size());

  for (unsigned int i = 0, n = chii.size(); i < n; ++i)
    xs[i] = tan(0.5*(chi0 - chii[i]));

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

  if (fRealAtm)
    MinuitFitErrors(chi0, t0, t0Err, rp, rpErr, chi2);
}

void 
ChiZeroRegression::MinuitFitFunc(int& /*npar*/, double* /*gin*/, double& f, double* par, int /*iflag*/)
{
  double& Chi0 = par[0];
  double& T0 = par[1];
  double& Rp = par[2];
  fdet::FTimeFitModel& timeFitModel = fdet::FTimeFitModel::GetInstance();

  double logLikeSum = 0.;
  for (int i = 0, n = chii.size(); i < n; ++i) {
    double ti_exp = timeFitModel.GetTimeAtAperture(T0, Rp, Chi0, chii[i], fsdpTheta, fEyeId, fTelId);
    logLikeSum += (-0.5*log(timeErrs[i]*timeErrs[i]) - 0.5*pow( (times[i] - ti_exp)/timeErrs[i] , 2));
  }
  f = -2.*logLikeSum;
}

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