testQuadraticFitter.cc

Go to the documentation of this file.

#include <utl/Histogram.h>
#include <utl/QuadraticFitter.h>
#include <tst/Verify.h>
#include <cppunit/extensions/HelperMacros.h>
#include <iostream>

using namespace std;
using namespace utl;
using namespace tst;


#define EQUAL(_x_, _y_) CPPUNIT_ASSERT(Verify<Equal>(_x_, _y_))
#define CLOSE(_x_, _y_) CPPUNIT_ASSERT(Verify<CloseTo>(_x_, _y_))
#define CLOSEEPS(_x_, _y_, _eps_) CPPUNIT_ASSERT(Verify<CloseTo>(_x_, _y_, _eps_))


class TestQuadraticFitter : public CppUnit::TestFixture {

  CPPUNIT_TEST_SUITE(TestQuadraticFitter);
  CPPUNIT_TEST(testHistogramExactFit);
  CPPUNIT_TEST(testHistogramSemiRandomFit);
  CPPUNIT_TEST(testFitRange);
  CPPUNIT_TEST_SUITE_END();

public:
  void
  setUp()
  {
    fA = -3;
    fB = 10000;
    fExtreme = 20;

    fAExp = 600;
    fDecay = 30;
    fExtreme2 = 50;
    fA2 = -1;
    fB2 = 400;
  }

  void tearDown() { }

  double quadraticFunction(const double x) const
  { return fA * Sqr(x - fExtreme) + fB; }

  double
  twoHumpFunction(const double x)
    const
  {
    const double par = fA2 * Sqr(x - fExtreme2) + fB2;
    const double dec = fAExp * exp(-x/fDecay);
    return dec + (par >= 0 ? par : 0);
  }

  void
  testHistogramExactFit()
  {
    const int n = 50;
    vector<double> binEdges;
    vector<double> counts;
    for (int i = 0; i < n; ++i) {
      binEdges.push_back(i);
      counts.push_back(quadraticFunction(i+0.5));
    }
    binEdges.push_back(n);

    VariableBinHistogramWrap<double, double> h(binEdges, counts);

    const double pos =
      MakeQuadraticFitter(h, FlatErrors()).GetExtremePosition();

    CLOSE(pos, fExtreme);

    double posErr;
    int ndof;
    const double pos2 =
      MakeQuadraticFitter(h, FlatErrors()).GetExtremePosition(posErr, ndof);

    CLOSE(pos2, fExtreme);
    EQUAL(ndof, 47);

    double chi2;
    const double pos3 =
      MakeQuadraticFitter(h, FlatErrors()).GetExtremePosition(posErr, chi2, ndof);

    CLOSE(pos3, fExtreme);
    CLOSE(chi2, 0.);
    EQUAL(ndof, 47);

    QuadraticFitData qf;
    MakeQuadraticFitter(h, FlatErrors()).GetFitData(qf);

    CLOSE(qf.GetExtremePosition(), fExtreme);
    CLOSE(qf.GetCoefficients()[0], fB);
    CLOSE(qf.GetCoefficients()[1], 0.);
    CLOSE(qf.GetCoefficients()[2], fA);

    CLOSE(qf(fExtreme), fB);
    CLOSE(qf(fExtreme + 3), quadraticFunction(fExtreme + 3));
    CLOSE(qf(fExtreme - 5), quadraticFunction(fExtreme - 5));
  }

  void
  testHistogramSemiRandomFit()
  {
    const int n = 50;
    vector<double> binEdges(n+1);
    vector<double> counts(n);
    for (int i = 0; i < n; ++i) {
      binEdges[i] = i;
      const double q = quadraticFunction(i+0.5);
      const double sq = q > 0 ? sqrt(2*q) : 0;
      const double qq = q + sq*sin(123456789.*i);
      counts[i] = qq;
    }
    binEdges[n] = n;

    VariableBinHistogramWrap<double, double> h(binEdges, counts);

    double posErr;
    int ndof;
    double chi2;
    const double pos =
      MakeQuadraticFitter(h).GetExtremePosition(posErr, chi2, ndof);

    cout << "\n pos=" << pos << " posErr=" << posErr << " chi2=" << chi2 << " "
            "ndof=" << ndof << endl;

    CLOSEEPS(pos, fExtreme, posErr);
    CLOSEEPS(posErr, 0., 1.);
    EQUAL(ndof, 47);

    const double omega = (chi2 - ndof) / sqrt(2.*ndof);

    CLOSEEPS(omega, 0., 0.5);
  }

  void
  testFitRange()
  {
    typedef unsigned short int ushort;
    const int n = 80;
    vector<char> binEdges(n+1);
    vector<ushort> counts(n);
    for (char i = 0; i < n; ++i) {
      binEdges[i] = i;
      counts[i] = ushort(twoHumpFunction(i));
    }
    binEdges[n] = char(n);

    VariableBinHistogramWrap<char, ushort> h(binEdges, counts);

    const int start = 30;
    const int stop = 70;

    double posErr;
    int ndof;
    double chi2;
    const double pos =
      MakeQuadraticFitter(h, start, stop).GetExtremePosition(posErr, chi2, ndof);

    cout << "\n pos=" << pos << " posErr=" << posErr << " chi2=" << chi2 << " "
            "ndof=" << ndof << endl;

    double posErr2;
    int ndof2;
    double chi22;
    const double pos2 =
      MakeQuadraticFitter(h, double(start), double(stop)).GetExtremePosition(posErr2, chi22, ndof2);

    CLOSE(pos, pos2);
    CLOSE(chi2, chi22);
    CLOSE(pos, pos2);
    EQUAL(ndof, ndof2);
    const double numericalSolution = 47.979666;
    CLOSEEPS(pos, numericalSolution + 0.5, posErr/10);

    cout << "\n pos=" << pos << " posErr=" << posErr << " chi2=" << chi2 << " "
            "ndof=" << ndof << endl;

    cout << "start=" << h.GetBinIndex(30.) << " stop=" << h.GetBinIndex(70.) << endl;
  }

private:
  double fA = 0;
  double fB = 0;
  double fExtreme = 0;

  double fAExp = 0;
  double fDecay = 0;
  double fExtreme2 = 0;
  double fA2 = 0;
  double fB2 = 0;
};


CPPUNIT_TEST_SUITE_REGISTRATION(TestQuadraticFitter);