GausFitCache.h

Go to the documentation of this file.

#ifndef GAUSFITCACHE_H
#define GAUSFITCACHE_H 1

// C/C++ standard libraries
#include <string>
#include <vector>

// ROOT libraries
#include "RtypesCore.h" // Double_t
#include "TF1.h"

namespace hit {

  class GausFitCache {
  public:
    GausFitCache(std::string new_name = "GausFitCache") : name(new_name) {}

    virtual ~GausFitCache();

    std::string GetName() const { return name; }

    virtual TF1* Get(size_t nFunc);

    virtual TF1* GetClone(size_t nGaus);

    virtual std::string FunctionName(size_t nFunc) const;

  protected:
    std::string name; 

    std::vector<TF1*> funcs;

    virtual TF1* CreateFunction(size_t nFunc) const;

  }; // class GausFitCache

  namespace details {

    template <typename T>
    inline T sqr(T v)
    {
      return v * v;
    }

    template <unsigned int NArg, typename FirstArg, typename... Args>
    struct TemplateArgumentHelper {
      using type = typename TemplateArgumentHelper<NArg - 1, Args...>::type;
    }; // struct TemplateArgumentHelper

    template <unsigned int NArg, typename... Args>
    struct TemplateArgument {
      using type = typename TemplateArgumentHelper<NArg, Args...>::type;
    }; // struct TemplateArgument

    template <unsigned int NFunc,
              Double_t Func(Double_t const*, Double_t const*),
              unsigned int NFuncParams>
    struct FuncSum {
      static Double_t eval(Double_t const*, Double_t const*);
      static constexpr unsigned int NParams = NFunc * NFuncParams;
    }; // struct FuncSum

    // partial specialization (declaration)
    template <Double_t Func(Double_t const*, Double_t const*), unsigned int NFuncParams>
    struct FuncSum<0U, Func, NFuncParams> {
      static Double_t eval(Double_t const*, Double_t const*);
      static constexpr unsigned int NParams = 0;
    }; // struct FuncSum<0>

    class CompiledGausFitCacheBaseStruct : public GausFitCache {
      /*
       * Since, as most of C++ metaprogramming code, this one is quite messy,
       * some explanations follow.
       *
       * The idea is to have a function sum of N base functions, and both N and
       * which base function can be specified as template parameters.
       *
       * The implementation of the sum of N base functions requires a
       * compile-time "for loop" that is implemented by recursion: a function
       * with a template parameter N=i calls the same with N=i-1. The loop is
       * interrupted by a template specialization for N=0 that does not call
       * N=-1. Unfortunately there is another template parameter that is carried
       * around, that is the base function. This means that the specialization
       * for N=0 is a partial specialization, since the function still needs
       * to be a template parameter. Partial specializations are only allowed
       * for classes (not for functions). Therefore we have here a class
       * with two template parameters and a static member eval() where the
       * total function is defined. Partial specializations of the full class
       * will implement the necessary N=0 exception to the recursion rule.
       *
       * Now, we want to put these in a vector of functions.
       * This is another loop, and since the function types are decided at
       * compile time, the loop must be done at compile time too.
       * Here it goes another class, moving around a vector.
       * Furthermore, now the actual type of function is a different one for
       * each element in the array: in this compile-time for-loop, the type
       * depends on the running parameter, and potentially on others.
       *
       *
       *
       *
       */
    public:
      using GausFitCache::GausFitCache;

      virtual TF1* GetClone(size_t nGaus);

      virtual unsigned int MaxGaussians() const { return funcs.size() - 1; }

      static Double_t gaus(Double_t const* x, Double_t const* params);

      template <unsigned int CutOff>
      static Double_t gaus_trunc(Double_t const* x, Double_t const* params);

      template <unsigned int NGaus>
      static Double_t ngaus(Double_t const* x, Double_t const* params)
      {
        return gaus(x, params) + ngaus<NGaus - 1>(x, params + 3);
      }

      template <unsigned int NGaus, unsigned int CutOff>
      static Double_t ngaus_trunc(Double_t const* x, Double_t const* params)
      {
        return FuncSum<NGaus, gaus_trunc<CutOff>, 3U>::eval(x, params);
      }

      template <unsigned int NGaus, unsigned int CutOff>
      using NGaussTruncClass = FuncSum<NGaus, gaus_trunc<CutOff>, 3U>;

    protected:
      template <unsigned int NFunc, template <unsigned int> class Func>
      struct InitializeFuncSumVector {
        static void fill(CompiledGausFitCacheBaseStruct& cache);
      }; // struct InitializeFuncSumVector

      // partial specialization: 0 of any function
      template <template <unsigned int> class Func>
      struct InitializeFuncSumVector<0U, Func> {
        static void fill(CompiledGausFitCacheBaseStruct& cache);
      }; // struct InitializeFuncSumVector<0, Func>

      template <unsigned int NGaus>
      void InitializeCompiledGausFitVector();

      template <unsigned int NGaus>
      void AppendFunction();

      void CannotCreateFunction [[noreturn]] (size_t nGaus) const;

    }; // class CompiledGausFitCacheBaseStruct

  } // namespace details

  template <unsigned int MaxGaus = 10>
  class CompiledGausFitCache : public details::CompiledGausFitCacheBaseStruct {
  public:
    CompiledGausFitCache(std::string new_name = "CompiledGausFitCache")
      : details::CompiledGausFitCacheBaseStruct(new_name)
    {
      InitializeCompiledGausFitVector<MaxGaus>();
    }

    virtual unsigned int MaxGaussians() const { return StoredMaxGaussians(); }

    constexpr unsigned int StoredMaxGaussians() const { return MaxGaus; }

  protected:
    virtual TF1* CreateFunction [[noreturn]] (size_t nGaus) const { CannotCreateFunction(nGaus); }

  }; // class CompiledGausFitCache

  template <unsigned int MaxGaus = 10, unsigned int CutOff = 5>
  class CompiledTruncatedGausFitCache : public details::CompiledGausFitCacheBaseStruct {
    template <unsigned int NGaus>
    using CutOffNGaussianClass = NGaussTruncClass<NGaus, CutOff>;

  public:
    CompiledTruncatedGausFitCache(std::string new_name = "CompiledTruncatedGausFitCache")
      : details::CompiledGausFitCacheBaseStruct(new_name)
    {
      InitializeFuncSumVector<MaxGaus, CutOffNGaussianClass>::fill(*this);
    }

    virtual unsigned int MaxGaussians() const { return StoredMaxGaussians(); }

    constexpr unsigned int StoredMaxGaussians() const { return MaxGaus; }

  protected:
    virtual TF1* CreateFunction [[noreturn]] (size_t nGaus) const { CannotCreateFunction(nGaus); }

  }; // class CompiledTruncatedGausFitCache

  //
  // template implementation
  //

  namespace details {

    // --- FuncSum -------------------------------------------------------------
    template <unsigned int NFunc,
              Double_t Func(Double_t const*, Double_t const*),
              unsigned int NFuncParams>
    constexpr unsigned int FuncSum<NFunc, Func, NFuncParams>::NParams;

    template <unsigned int NFunc,
              Double_t Func(Double_t const*, Double_t const*),
              unsigned int NFuncParams>
    Double_t FuncSum<NFunc, Func, NFuncParams>::eval(Double_t const* x, Double_t const* params)
    {
      return Func(x, params + NFuncParams * (NFunc - 1)) // use the last parameters
             + FuncSum<NFunc - 1, Func, NFuncParams>::eval(x, params);
    } // CompiledGausFitCacheBaseStruct::FuncSum<NFunc, Func>::eval()

    // partial specialization: 0 of any function
    template <Double_t Func(Double_t const*, Double_t const*), unsigned int NFuncParams>
    Double_t FuncSum<0U, Func, NFuncParams>::eval(Double_t const*, Double_t const*)
    {
      return 0.;
    }

    // --- CompiledGausFitCacheBaseStruct --------------------------------------

    template <unsigned int NGaus>
    void CompiledGausFitCacheBaseStruct::InitializeCompiledGausFitVector()
    {
      if (NGaus > 0) InitializeCompiledGausFitVector<NGaus - 1>();
      AppendFunction<NGaus>();
    } // CompiledGausFitCacheBaseStruct::InitializeCompiledGausFitVector()

    template <>
    inline void CompiledGausFitCacheBaseStruct::InitializeCompiledGausFitVector<0>()
    {
      AppendFunction<0>();
    }

    template <unsigned int NGaus>
    void CompiledGausFitCacheBaseStruct::AppendFunction()
    {
      // create a function in the ficticious range [ 0, 1 ]:
      funcs.push_back(new TF1(FunctionName(NGaus).c_str(), &ngaus<NGaus>, 0., 1., 3 * NGaus));
    } // CompiledGausFitCacheBaseStruct::AppendFunction()

    template <unsigned int CutOff>
    Double_t CompiledGausFitCacheBaseStruct::gaus_trunc(Double_t const* x, Double_t const* params)
    {
      const Double_t z = (x[0] - params[1]) / params[2];
      return ((z > -((Double_t)CutOff)) && (z < (Double_t)CutOff)) ?
               params[0] * std::exp(-0.5 * sqr(z)) :
               0.;
    } // CompiledGausFitCacheBaseStruct::gaus_trunc()

    template <>
    inline Double_t CompiledGausFitCacheBaseStruct::ngaus<0>(Double_t const* /* x */,
                                                             Double_t const* /* params */)
    {
      return 0.;
    }

    // --- CompiledGausFitCacheBaseStruct::InitializeFuncSumVector -------------

    template <unsigned int NFunc, template <unsigned int> class Func>
    void CompiledGausFitCacheBaseStruct::InitializeFuncSumVector<NFunc, Func>::fill(
      CompiledGausFitCacheBaseStruct& cache)
    {
      // first fill the lower functions
      InitializeFuncSumVector<NFunc - 1, Func>::fill(cache);
      // then add one
      cache.funcs.push_back(new TF1(
        cache.FunctionName(NFunc).c_str(), Func<NFunc>::eval, 0., 1., Func<NFunc>::NParams));
    } // InitializeFuncSumVector<NFunc, Func>::fill()

    template <template <unsigned int> class Func>
    void CompiledGausFitCacheBaseStruct::InitializeFuncSumVector<0U, Func>::fill(
      CompiledGausFitCacheBaseStruct& cache)
    {
      cache.funcs.push_back(new TF1(cache.FunctionName(0).c_str(), Func<0U>::eval, 0., 1., 0));
    } // InitializeFuncSumVector<0, Func>::fill()

    // -------------------------------------------------------------------------

  } // namespace details

  // ---------------------------------------------------------------------------

} // namespace hit

#endif // GAUSFITCACHE_H