ROptFitter.cc

Go to the documentation of this file.

#include "ROptFitter.h"

#include <Minuit2/FCNBase.h>
#include <Minuit2/MnMinimize.h>
#include <Minuit2/MnHesse.h>
#include <Minuit2/FunctionMinimum.h>
#include <Minuit2/MnUserParameters.h>

#include <utl/AugerUnits.h>
#include <utl/Math.h>

#include <cmath>
#include <vector>

using namespace LDFFinderKG;


bool
ROptFitter::GetResult(double& rOpt)
  const
{
  ROOT::Minuit2::MnUserParameters pars;
  pars.Add("rOpt", 1400*utl::meter, 50*utl::meter);
  pars.SetLowerLimit("rOpt", 0);

  ROOT::Minuit2::MnMinimize m(*this, pars, 0);

  // TODO make minuit verbose

  ROOT::Minuit2::FunctionMinimum fmin = m();

  ROOT::Minuit2::MnHesse hesse;
  hesse(m.Fcnbase(), fmin);

  if (!fmin.HasValidParameters())
    return false;

  rOpt = fmin.UserParameters().Params()[0];
  return true;
}


double
ROptFitter::operator()(const std::vector<double>& pars)
  const
{
  const double& ropt = pars[0];

  const int n = fShowerSize.size();
  double sum = 0;
  double sum2 = 0;
  for (int i = 0; i < n; ++i) {
    fLDFParameters[0] = fBeta[i];
    const double s = fShowerSize[i] * fLDFFunction(ropt, fLDFParameters);
    sum += s;
    sum2 += utl::Sqr(s);
  }
  const double invN = 1./n;
  const double avg = sum * invN;
  const double var2 = sum2 * invN - utl::Sqr(avg);

  return var2 > 0 ? std::sqrt(var2)/avg : 0;
}