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 "Rtypes.h" // Double_t
#include "TF1.h" // Double_t

// Framework libraries
#include "canvas/Utilities/Exception.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