rogly_fitting_function.h

#ifndef ROGLY_FITTING_FUNCTION_H
#define ROGLY_FITTING_FUNCTION_H
 
#include "ipde_fit_parameter.h"
#include <TH1.h>
#include <vector>
#include "TRandom.h"
 
namespace KVImpactParameters {
   template<unsigned int PolyDegree>
   class rogly_fitting_function {
      ipde_fit_parameter theta_p;
      ipde_fit_parameter kmax;
      std::vector<ipde_fit_parameter> poly_param;
 
   public:
      rogly_fitting_function()
         : poly_param(PolyDegree)
      {}
      rogly_fitting_function(const rogly_fitting_function& prev_fit)
         : theta_p(prev_fit.theta_p), kmax(prev_fit.kmax), poly_param(PolyDegree)
      {
         poly_param = prev_fit.poly_param;
      }
      rogly_fitting_function(std::initializer_list<double> params)
         : poly_param(PolyDegree)
      {
         int i = 0;
         for (auto param : params) {
            switch (i) {
               case 0:
                  theta_p.value = param;
                  break;
               case 1:
                  kmax.value = param;
                  break;
               default:
                  poly_param[i - 2].value = param;
            }
            ++i;
         }
      }
      constexpr int npar() const
      {
         return PolyDegree + 2;
      }
 
      double k_cb(double cb) const
      {
         double kcb = 0;
         for (int j = 1; j <= PolyDegree; ++j) kcb += poly_param[j - 1].value * pow(cb, j);
         return kmax.value * TMath::Exp(-kcb);
      }
      void fill_params_from_array(double* p)
      {
         theta_p.value = p[0];
         kmax.value = p[1];
         for (int j = 1; j <= PolyDegree; ++j) poly_param[j - 1].value = p[j + 1];
      }
      void fill_array_from_params(double* p) const
      {
         p[0] = theta_p.value;
         p[1] = kmax.value;
         for (int j = 1; j <= PolyDegree; ++j) p[j + 1] = poly_param[j - 1].value ;
      }
      void set_par_names(TF1& f) const
      {
         f.SetParName(0, "#theta");
         f.SetParName(1, "k_{max}");
         for (int j = 1; j <= PolyDegree; ++j) f.SetParName(j + 1, Form("a_{%d}", j));
      }
      void set_initial_parameters(TH1* h, TF1& f)
      {
         f.SetParameter(0, theta_p.value = 0.5);
         f.SetParLimits(0, 0.1, 5);
         double xmax = h->GetXaxis()->GetBinUpEdge(h->GetNbinsX());
         double xmin = xmax / 7.;
         f.SetParameter(1, kmax.value = (xmin * 2) / theta_p.value);
         f.SetParLimits(1, xmin / theta_p.value, xmax / theta_p.value);
         for (int j = 1; j <= PolyDegree; ++j) {
            f.SetParameter(j + 1, gRandom->Uniform(-1, 1));
            f.SetParLimits(j + 1, -5, 5);
         }
      }
      void print_fit_params() const
      {
         std::cout << "Theta = ";
         theta_p.print();
         std::cout << "\n";
         std::cout << "kmax = ";
         kmax.print();
         std::cout << "(Xmax = " << theta_p.value* kmax.value << ")";
         std::cout << "\n";
         for (int j = 1; j <= PolyDegree; ++j) {
            std::cout << "a" << j << " = ";
            poly_param[j - 1].print();
            std::cout << "  ";
         }
         std::cout << "\n";
      }
      void backup_params()
      {
         theta_p.backup();
         kmax.backup();
         std::for_each(poly_param.begin(), poly_param.end(), [](ipde_fit_parameter & p) {
            p.backup();
         });
      }
      void restore_params()
      {
         theta_p.restore();
         kmax.restore();
         std::for_each(poly_param.begin(), poly_param.end(), [](ipde_fit_parameter & p) {
            p.restore();
         });
      }
      void normalise_shape_function()
      {
         kmax.value = 1;
      }
      double meanX(double cb) const
      {
         return theta_p.value * k_cb(cb);
      }
      double redVar(double) const
      {
         return theta_p.value;
      }
 
      ClassDef(rogly_fitting_function, 0) //fit using exponential polynomial
   };
}
 
#endif // ROGLY_FITTING_FUNCTION_H