RoptFit.cc

Go to the documentation of this file.

#include <TMinuit.h>
#include <utl/AugerUnits.h>
#include "RoptFit.h"

using namespace LDFFinderOG;
using namespace std;


LDFFunctionType RoptFit::fgLDFType;
int RoptFit::fgN;
const double* RoptFit::fgS1000;
const double* RoptFit::fgBeta;


pair<double, double>
RoptFit::Fit(const LDFFunctionType type,
             const int n, const double* const s1000, const double* const beta,
             const bool verbose)
{
  fgLDFType = type;
  fgN = n;
  fgS1000 = s1000;
  fgBeta = beta;

  TMinuit minuit(1);
  int errFlag = 0;
  double argList[10];
  argList[0] = -1;
  if (!verbose) {
    minuit.mnexcm("SET PRINTOUT", argList, 1, errFlag);
    minuit.mnexcm("SET NOWARNINGS", argList, 0, errFlag);
    minuit.SetPrintLevel(-1);
  }

  minuit.SetFCN(RoptFit::FitFunction);
  argList[0] = 1;
  minuit.mnexcm("SET ERRORDEF", argList, 1, errFlag);

  minuit.mnparm(0, "ropt", 1400*meter, 50*meter, 0, 0, errFlag);

  argList[0] = 500;
  minuit.mnexcm("MINIMIZE", argList, 1, errFlag);

  // Get results
  pair<double, double> ropt;
  minuit.GetParameter(0, ropt.first, ropt.second);

  return ropt;
}


void
RoptFit::FitFunction(int& /*nPar*/, double* const /*grad*/, double& value,
                     double* const par, const int /*flag*/)
{
  const double& ropt = par[0];

  double sum = 0;
  double sum2 = 0;
  for (int i = 0; i < fgN; ++i) {
    const double s = fgS1000[i] * LDFFunction(fgLDFType, ropt, fgBeta[i], 0);
    sum += s;
    sum2 += Sqr(s);
  }
  const double invN = 1./fgN;
  const double avg = sum * invN;
  const double var2 = sum2 * invN - Sqr(avg);

  value = (var2 > 0 ? sqrt(var2) : 0) / avg;
}


// Configure (x)emacs for this file ...
// Local Variables:
// mode: c++
// End: