GausFitCache_test.cc File Reference

Test for classes in GausFitCache,h. More...

#include <cassert>
#include <cmath>
#include <array>
#include <limits>
#include "boost/test/unit_test.hpp"
#include "cetlib/pow.h"
#include "TH1D.h"
#include "TFitResultPtr.h"
#include "TFitResult.h"
#include "larreco/RecoAlg/GausFitCache.h"

Go to the source code of this file.

Classes
struct	RootGausFuncWrapper

Macros
#define	BOOST_TEST_MODULE   ( HitAnaAlg_test )

Typedefs
using	tolerance_t = decltype(0.001%tolerance())

Functions
	BOOST_AUTO_TEST_CASE (TestGaussianTest)
	BOOST_AUTO_TEST_CASE (GaussianTest)
	BOOST_AUTO_TEST_CASE (GaussianTrunc5Test)
	BOOST_AUTO_TEST_CASE (GaussianTrunc4Test)
Double_t	multi_gaus (Double_t x, const unsigned int nGaus, Double_t const *params)
	Expect for each Gaussian ROOT-like parameters: amplitude, mean, sigma. More...
std::vector< Double_t >	SortGaussianResults (TF1 const *pFunc, Double_t const *Params)
void	ThreeGaussianFitTest (hit::GausFitCache &GausCache, tolerance_t tol)
	BOOST_AUTO_TEST_CASE (ThreeGaussianTest)
	BOOST_AUTO_TEST_CASE (RunTimeThreeGaussianFitTest)
	BOOST_AUTO_TEST_CASE (CompiledThreeGaussianFitTest)
	BOOST_AUTO_TEST_CASE (CompiledTruncated5ThreeGaussianFitTest)
	BOOST_AUTO_TEST_CASE (CompiledTruncated4ThreeGaussianFitTest)
	BOOST_AUTO_TEST_CASE (CompiledTruncated3ThreeGaussianFitTest)

Detailed Description

Test for classes in GausFitCache,h.

Author

Gianluca Petrillo (petri.nosp@m.llo@.nosp@m.fnal..nosp@m.gov)

Date

March 27th, 2015

See also

GausFitCache,h

Definition in file GausFitCache_test.cc.

Macro Definition Documentation

#define BOOST_TEST_MODULE   ( HitAnaAlg_test )

Definition at line 19 of file GausFitCache_test.cc.

Typedef Documentation

using tolerance_t = decltype(0.001% tolerance())

Definition at line 36 of file GausFitCache_test.cc.

Function Documentation

BOOST_AUTO_TEST_CASE

(

TestGaussianTest

)

Definition at line 52 of file GausFitCache_test.cc.

{
  // static function to be tested:
  const auto gaus_func = ::gaus;

  const double exp1sigma = std::exp(-0.5*square(1.));
  const double exp2sigma = std::exp(-0.5*square(2.));

  // here there should be some more believable test...
  const double amplitude = 16.;
  for (double sigma = 0.5; sigma < 2.25; ++sigma) {
    for (double mean = -5; mean < 5.5; ++mean) {

      BOOST_TEST(gaus_func(-6.5 * sigma + mean, mean, sigma, amplitude) > 0.);
      BOOST_TEST(gaus_func(-5.5 * sigma + mean, mean, sigma, amplitude) > 0.);
      BOOST_TEST(gaus_func(-4.5 * sigma + mean, mean, sigma, amplitude) > 0.);

      auto const tol = 0.001% tolerance();
      BOOST_TEST(gaus_func(-2.0 * sigma + mean, mean, sigma, amplitude) == amplitude * exp2sigma, tol);
      BOOST_TEST(gaus_func(-1.0 * sigma + mean, mean, sigma, amplitude) == amplitude * exp1sigma, tol);
      BOOST_TEST(gaus_func( 0.0 * sigma + mean, mean, sigma, amplitude) == amplitude, tol);
      BOOST_TEST(gaus_func( 1.0 * sigma + mean, mean, sigma, amplitude) == amplitude * exp1sigma, tol);
      BOOST_TEST(gaus_func( 2.0 * sigma + mean, mean, sigma, amplitude) == amplitude * exp2sigma, tol);

      BOOST_TEST(gaus_func( 4.5 * sigma + mean, mean, sigma, amplitude) > 0.);
      BOOST_TEST(gaus_func( 5.5 * sigma + mean, mean, sigma, amplitude) > 0.);
      BOOST_TEST(gaus_func( 6.5 * sigma + mean, mean, sigma, amplitude) > 0.);

    } // for mean
  } // for sigma

} // BOOST_AUTO_TEST_CASE(TestGaussianTest)

BOOST_AUTO_TEST_CASE

(

GaussianTest

)

Definition at line 107 of file GausFitCache_test.cc.

{
  // static function to be tested:
  const auto gaus_func = RootGausFuncWrapper
    (hit::details::CompiledGausFitCacheBaseStruct::gaus);

  const double exp1sigma = std::exp(-0.5*square(1.));
  const double exp2sigma = std::exp(-0.5*square(2.));

  const double amplitude = 16.;
  for (double sigma = 0.5; sigma < 2.25; ++sigma) {
    for (double mean = -5; mean < 5.5; ++mean) {

      BOOST_TEST(gaus_func(-6.5 * sigma + mean, mean, sigma, amplitude) > 0.);
      BOOST_TEST(gaus_func(-5.5 * sigma + mean, mean, sigma, amplitude) > 0.);
      BOOST_TEST(gaus_func(-4.5 * sigma + mean, mean, sigma, amplitude) > 0.);

      auto const tol = 0.001% tolerance();
      BOOST_TEST(gaus_func(-2.0 * sigma + mean, mean, sigma, amplitude) == amplitude * exp2sigma, tol);
      BOOST_TEST(gaus_func(-1.0 * sigma + mean, mean, sigma, amplitude) == amplitude * exp1sigma, tol);
      BOOST_TEST(gaus_func( 0.0 * sigma + mean, mean, sigma, amplitude) == amplitude, tol);
      BOOST_TEST(gaus_func(+1.0 * sigma + mean, mean, sigma, amplitude) == amplitude * exp1sigma, tol);
      BOOST_TEST(gaus_func(+2.0 * sigma + mean, mean, sigma, amplitude) == amplitude * exp2sigma, tol);

      BOOST_TEST(gaus_func(+4.5 * sigma + mean, mean, sigma, amplitude) > 0.);
      BOOST_TEST(gaus_func(+5.5 * sigma + mean, mean, sigma, amplitude) > 0.);
      BOOST_TEST(gaus_func(+6.5 * sigma + mean, mean, sigma, amplitude) > 0.);

    } // for mean
  } // for sigma

} // BOOST_AUTO_TEST_CASE(GaussianTest)

BOOST_AUTO_TEST_CASE

(

GaussianTrunc5Test

)

Definition at line 141 of file GausFitCache_test.cc.

{
  // static function to be tested:
  const auto gaus_func = RootGausFuncWrapper
    (hit::details::CompiledGausFitCacheBaseStruct::gaus_trunc<5>);

  const double exp1sigma = std::exp(-0.5*square(1.));
  const double exp2sigma = std::exp(-0.5*square(2.));

  const double amplitude = 16.;
  for (double sigma = 0.5; sigma < 2.25; ++sigma) {
    for (double mean = -5; mean < 5.5; ++mean) {

      BOOST_TEST(gaus_func(-6.5 * sigma + mean, mean, sigma, amplitude) == 0.);
      BOOST_TEST(gaus_func(-5.5 * sigma + mean, mean, sigma, amplitude) == 0.);
      BOOST_TEST(gaus_func(-4.5 * sigma + mean, mean, sigma, amplitude) > 0.);

      auto const tol = 0.001% tolerance();
      BOOST_TEST(gaus_func(-2.0 * sigma + mean, mean, sigma, amplitude) == amplitude * exp2sigma, tol);
      BOOST_TEST(gaus_func(-1.0 * sigma + mean, mean, sigma, amplitude) == amplitude * exp1sigma, tol);
      BOOST_TEST(gaus_func( 0.0 * sigma + mean, mean, sigma, amplitude) == amplitude, tol);
      BOOST_TEST(gaus_func( 1.0 * sigma + mean, mean, sigma, amplitude) == amplitude * exp1sigma, tol);
      BOOST_TEST(gaus_func( 2.0 * sigma + mean, mean, sigma, amplitude) == amplitude * exp2sigma, tol);

      BOOST_TEST(gaus_func( 4.5 * sigma + mean, mean, sigma, amplitude) > 0.);
      BOOST_TEST(gaus_func( 5.5 * sigma + mean, mean, sigma, amplitude) == 0.);
      BOOST_TEST(gaus_func( 6.5 * sigma + mean, mean, sigma, amplitude) == 0.);

    } // for mean
  } // for sigma

} // BOOST_AUTO_TEST_CASE(GaussianTrunc5Test)

BOOST_AUTO_TEST_CASE

(

GaussianTrunc4Test

)

Definition at line 175 of file GausFitCache_test.cc.

{
  // static function to be tested:
  const auto gaus_func = RootGausFuncWrapper
    (hit::details::CompiledGausFitCacheBaseStruct::gaus_trunc<4>);

  const double exp1sigma = std::exp(-0.5*square(1.));
  const double exp2sigma = std::exp(-0.5*square(2.));

  const double amplitude = 16.;
  for (double sigma = 0.5; sigma < 2.25; ++sigma) {
    for (double mean = -5; mean < 5.5; ++mean) {

      BOOST_TEST(gaus_func(-6.5 * sigma + mean, mean, sigma, amplitude) == 0.);
      BOOST_TEST(gaus_func(-5.5 * sigma + mean, mean, sigma, amplitude) == 0.);
      BOOST_TEST(gaus_func(-4.5 * sigma + mean, mean, sigma, amplitude) == 0.);

      auto const tol = 0.001% tolerance();
      BOOST_TEST(gaus_func(-2.0 * sigma + mean, mean, sigma, amplitude) == amplitude * exp2sigma, tol);
      BOOST_TEST(gaus_func(-1.0 * sigma + mean, mean, sigma, amplitude) == amplitude * exp1sigma, tol);
      BOOST_TEST(gaus_func( 0.0 * sigma + mean, mean, sigma, amplitude) == amplitude, tol);
      BOOST_TEST(gaus_func( 1.0 * sigma + mean, mean, sigma, amplitude) == amplitude * exp1sigma, tol);
      BOOST_TEST(gaus_func( 2.0 * sigma + mean, mean, sigma, amplitude) == amplitude * exp2sigma, tol);

      BOOST_TEST(gaus_func( 4.5 * sigma + mean, mean, sigma, amplitude) == 0.);
      BOOST_TEST(gaus_func( 5.5 * sigma + mean, mean, sigma, amplitude) == 0.);
      BOOST_TEST(gaus_func( 6.5 * sigma + mean, mean, sigma, amplitude) == 0.);

    } // for mean
  } // for sigma

} // BOOST_AUTO_TEST_CASE(GaussianTrunc4Test)

BOOST_AUTO_TEST_CASE

(

ThreeGaussianTest

)

Definition at line 364 of file GausFitCache_test.cc.

{

  const Double_t Params[] = {
    4.0, -2.0, 1.5,
    5.0,  0.1, 0.3,
    2.0,  1.0, 0.5
  }; // Params[]

  // rounding errors on x on the look do not really matter
  for (double x = -10.; x <= +10.; x += 0.1) {
    const Double_t expected = multi_gaus(x, 3, Params);

    const Double_t tested
      = hit::details::CompiledGausFitCacheBaseStruct::ngaus<3>(&x, Params);

    BOOST_TEST(tested == expected, 0.001% tolerance());
  } // for x

} // BOOST_AUTO_TEST_CASE(ThreeGaussianTest)

BOOST_AUTO_TEST_CASE

(

RunTimeThreeGaussianFitTest

)

Definition at line 387 of file GausFitCache_test.cc.

{
  // max 20 Gaussians
  hit::GausFitCache GausCache("RunTimeGaussians");
  ThreeGaussianFitTest(GausCache, 0.02% tolerance());
} // BOOST_AUTO_TEST_CASE(RunTimeThreeGaussianFitTest)

BOOST_AUTO_TEST_CASE

(

CompiledThreeGaussianFitTest

)

Definition at line 396 of file GausFitCache_test.cc.

{
   // max 20 Gaussians
  hit::CompiledGausFitCache<20> GausCache("CompiledGaussians");
  ThreeGaussianFitTest(GausCache, 0.02% tolerance());
} // BOOST_AUTO_TEST_CASE(CompiledThreeGaussianFitTest)

BOOST_AUTO_TEST_CASE

(

CompiledTruncated5ThreeGaussianFitTest

)

Definition at line 406 of file GausFitCache_test.cc.

{
   // max 20 Gaussians; truncate at 5 sigma
  hit::CompiledTruncatedGausFitCache<20, 5> GausCache
    ("CompiledTruncated5Gaussians");
  ThreeGaussianFitTest(GausCache, 0.1% tolerance()); // 0.1% tolerance
} // BOOST_AUTO_TEST_CASE(CompiledTruncated5ThreeGaussianFitTest)

BOOST_AUTO_TEST_CASE

(

CompiledTruncated4ThreeGaussianFitTest

)

Definition at line 417 of file GausFitCache_test.cc.

{
   // max 20 Gaussians; truncate at 4 sigma
  hit::CompiledTruncatedGausFitCache<20, 4> GausCache
    ("CompiledTruncated4Gaussians");
  ThreeGaussianFitTest(GausCache, 0.1% tolerance());
} // BOOST_AUTO_TEST_CASE(CompiledTruncated4ThreeGaussianFitTest)

BOOST_AUTO_TEST_CASE

(

CompiledTruncated3ThreeGaussianFitTest

)

Definition at line 427 of file GausFitCache_test.cc.

{
   // max 20 Gaussians; truncate at 3 sigma
  hit::CompiledTruncatedGausFitCache<20, 3> GausCache
    ("CompiledTruncated3Gaussians");
  ThreeGaussianFitTest(GausCache, 10.% tolerance()); // (seriously, it's that bad)
} // BOOST_AUTO_TEST_CASE(CompiledTruncated3ThreeGaussianFitTest)

Double_t multi_gaus	(	Double_t	x,
		const unsigned int	nGaus,
		Double_t const *	params
	)

Expect for each Gaussian ROOT-like parameters: amplitude, mean, sigma.

Definition at line 214 of file GausFitCache_test.cc.

{

  Double_t res = 0.;
  for (unsigned int iGaus = 0; iGaus < nGaus; ++iGaus) {
    // our gaus() takes the parameters as: mean, sigma, amplitude
    res += gaus(x, params[1], params[2], params[0]);
    params += 3;
  } // for
  return res;
} // multi_gaus()

std::vector<Double_t> SortGaussianResults	(	TF1 const *	pFunc,
		Double_t const *	Params
	)

Definition at line 229 of file GausFitCache_test.cc.

{
  assert(pFunc->GetNpar() % 3 == 0); // Sorting non-Gaussian function!

  const size_t nGaus = pFunc->GetNpar() / 3;

  std::vector<size_t> BestMatch(nGaus, nGaus); // initialize with invalid number
  for (size_t iGaus = 0; iGaus < nGaus; ++iGaus) {
    double best_match_quality = std::numeric_limits<double>::max();
    for (size_t iFitGaus = 0; iFitGaus < nGaus; ++iFitGaus) {
      // so far, we compare all the parameters with the same weight
      const double match_quality
        = square(Params[iGaus*3 + 0] - pFunc->GetParameter(iFitGaus*3 + 0))
        + square(Params[iGaus*3 + 1] - pFunc->GetParameter(iFitGaus*3 + 1))
        + square(Params[iGaus*3 + 2] - pFunc->GetParameter(iFitGaus*3 + 2))
        ;
      if (match_quality < best_match_quality) {
        best_match_quality = match_quality;
        BestMatch[iGaus] = iFitGaus;
      } // if this is the best so far
    } // for fit

    assert(BestMatch[iGaus] < nGaus); // No reasonable solution matched!

    for (size_t i = 0; i < iGaus; ++i) {
      // Same fit Gaussian solution matches multiple solutions
      assert(BestMatch[i] != BestMatch[iGaus]);
    } // for i

  } // for solution

  std::vector<Double_t> sorted_params(nGaus * 3);
  for (size_t iGaus = 0; iGaus < nGaus; ++iGaus) {
    const size_t iBestFit = BestMatch[iGaus];
    for (size_t i = 0; i < 3; ++i)
      sorted_params[iGaus*3 + i] = pFunc->GetParameter(iBestFit*3 + i);
  }
  return sorted_params;
} // SortGaussianResults()

void ThreeGaussianFitTest	(	hit::GausFitCache &	GausCache,
		tolerance_t	tol
	)

Definition at line 271 of file GausFitCache_test.cc.

                                           {
  std::string name = GausCache.GetName();

  constexpr size_t nGaus = 3; // we test three gaussians

  // sorted by increasing mean
  const Double_t Params[nGaus * 3] = { // amplitude, mean, sigma
    4.0, -2.0, 1.5,
    5.0,  0.1, 0.3,
    2.0,  1.0, 0.5
  }; // Params[]

  // fit test
  // - fill the input histogram
  TH1D Hist
    ("H3Gaus", ("Three-Gaussian test - " + name).c_str(), 200, -10., 10.);
  for (Int_t iBin = 1; iBin <= Hist.GetNbinsX(); ++iBin)
    Hist.SetBinContent(iBin, multi_gaus(Hist.GetBinCenter(iBin), nGaus, Params));

#ifdef GAUSFITCACHE_TEST_DEBUG
  Hist.Draw();
#endif // GAUSFITCACHE_TEST_DEBUG

  // - get the 3-Gaussian fitting function
  TF1* pFunc = GausCache.Get(nGaus);

  // - prepare the function with reasonable initial values
  const double MeanOffset = -((double(nGaus) - 1.) / 2.);
  const double MeanStep = 1.;
  for (size_t iGaus = 0; iGaus < nGaus; ++iGaus) {
    size_t BaseIndex = iGaus * 3;
    pFunc->SetParameter(BaseIndex + 0, 1.); // amplitude
    pFunc->SetParameter(BaseIndex + 1, MeanOffset + iGaus * MeanStep); // mean
    pFunc->SetParameter(BaseIndex + 2, 1.); // sigma
  } // for

  // function range is ignored in the fit, but used when drawing
  pFunc->SetRange(Hist.GetBinCenter(1), Hist.GetBinCenter(Hist.GetNbinsX()));

#ifdef GAUSFITCACHE_TEST_DEBUG
  TF1* pFCopy = pFunc->DrawCopy("L SAME");
  pFCopy->SetLineStyle(kDashed);
  pFCopy->SetLineColor(kCyan);
#endif // GAUSFITCACHE_TEST_DEBUG

  // - fit
  TFitResultPtr fit_res = Hist.Fit(pFunc, "WQ0NS");
  BOOST_TEST(int(fit_res) == 0);

#ifdef GAUSFITCACHE_TEST_DEBUG
  pFCopy = pFunc->DrawCopy("L SAME");
  pFCopy->SetLineWidth(1);
  pFCopy->SetLineColor(kRed);

  gPad->SaveAs(("3GausTest_" + name + ".eps").c_str());
#endif // GAUSFITCACHE_TEST_DEBUG

  // - sort the results according to the expected ones
  std::vector<Double_t> results = SortGaussianResults(pFunc, Params);


  // - check them
  for (size_t iGaus = 0; iGaus < nGaus; ++iGaus) {
    const size_t iParam = iGaus * 3;

    // - check amplitude
    BOOST_TEST(results[iParam + 0] == Params[iParam + 0], tol);
    // - check mean
    BOOST_TEST(results[iParam + 1] == Params[iParam + 1], tol);
    // - check sigma
    BOOST_TEST(results[iParam + 2] == Params[iParam + 2], tol);

  } // for iGaus

  BOOST_WARN(pFunc->GetChisquare() < pFunc->GetNDF() / 2.);

#ifdef GAUSFITCACHE_TEST_DEBUG
  pFunc->Print();
  // always print the fit result
#else // !GAUSFITCACHE_TEST_DEBUG
  // in case of trouble, print the fit result
  if ((int(fit_res) != 0) || (pFunc->GetChisquare() >= pFunc->GetNDF() / 2.))
#endif // GAUSFITCACHE_TEST_DEBUG
  {
    fit_res->Print();
  }

} // ThreeGaussianFitTest()