AttenuationFit.h

Go to the documentation of this file.

// \author Darko Veberic
// \date 2014

#ifndef _cic_AttenuationFit_h_
#define _cic_AttenuationFit_h_

#include <Minuit2/FCNBase.h>
#include <Minuit2/MnMinimize.h>
#include <Minuit2/MnHesse.h>
#include <Minuit2/FunctionMinimum.h>
#include <Minuit2/MnUserParameters.h>


namespace cic {

  namespace AttenuationFit {

    typedef std::function<double(const double, const std::vector<double>&)> Attenuation;

    typedef std::function<double(double, double)> SFunction;


    template<class CIC>
    class CICFitter : public ROOT::Minuit2::FCNBase {
    public:
      explicit
      CICFitter(CIC& cic, const Attenuation& att, unsigned int ndof = 1)
        : fCIC(&cic), fAttenuation(att), fNdof(ndof) { }

      std::vector<double>
      Minimize(const std::vector<std::string>& name,
               const std::vector<double>& init,
               const std::vector<double>& sigma)
      {
        if (name.size() != init.size() ||
            name.size() != sigma.size())
          throw std::runtime_error("size mismatch!");

        ROOT::Minuit2::MnUserParameters pars;
        for (unsigned int i = 0, n = name.size(); i < n; ++i) {
          pars.Add(name[i], init[i], sigma[i]);
          if (!sigma[i])
            pars.Fix(name[i]);
        }

        ROOT::Minuit2::MnMinimize m(*this, pars, 0);
        ROOT::Minuit2::FunctionMinimum fmin = m();
        ROOT::Minuit2::MnHesse hesse;
        hesse(m.Fcnbase(), fmin);

        if (!fmin.HasValidParameters())
          return std::vector<double>();
        std::vector<double> result = fmin.UserParameters().Params();
        result.push_back(fmin.Fval());
        result.push_back(fNdof);
        return result;
      }

      double Up() const override { return 1; }

      virtual double operator()(const std::vector<double>& pars) const override = 0;

      double GetAttenuation(const double sin2theta, const std::vector<double>& pars) const
      { return fAttenuation(sin2theta, pars); }

    private:
      CICFitter(const CICFitter&) = delete;
      CICFitter& operator=(const CICFitter&) = delete;

    protected:
      CIC* fCIC = nullptr;
      Attenuation fAttenuation;
      unsigned int fNdof = 1;
    };


    template<class CIC>
    class KS : public CICFitter<CIC> {
    public:
      explicit
      KS(CIC& cic, const Attenuation& att, const double lnSCut)
        : CICFitter<CIC>(cic, att), fLnSCut(lnSCut) { }

      virtual
      double
      operator()(const std::vector<double>& pars)
        const
        override
      {
        const auto att = [this, &pars](const double sin2theta) {
          return this->GetAttenuation(sin2theta, pars);
        };
        return this->fCIC->GetKS(att, fLnSCut);
      }

    private:
      double fLnSCut = 0;
    };


    template<class CIC>
    class AndersonDarling : public CICFitter<CIC> {
    public:
      explicit
      AndersonDarling(CIC& cic, const Attenuation& att, const double lnSCut)
        : CICFitter<CIC>(cic, att), fLnSCut(lnSCut) { }

      virtual
      double
      operator()(const std::vector<double>& pars)
        const
        override
      {
        const auto att = [this, &pars](const double sin2theta) {
          return this->GetAttenuation(sin2theta, pars);
        };
        return this->fCIC->GetAndersonDarling(att, fLnSCut);
      }

    private:
      double fLnSCut = 0;
    };


    template<class CIC>
    class Chi2 : public CICFitter<CIC> {
    public:
      explicit
      Chi2(CIC& cic, const Attenuation& att, const double lnSCut)
        : CICFitter<CIC>(cic, att), fLnSCut(lnSCut) { }

      virtual
      double
      operator()(const std::vector<double>& pars)
        const
        override
      {
        const auto att = [this, &pars](const double sin2theta) {
          return this->GetAttenuation(sin2theta, pars);
        };
        return this->fCIC->GetChi2Ndof(att, fLnSCut);
      }

    private:
      double fLnSCut = 0;
    };


    template<class CIC>
    class PoissonChi2 : public CICFitter<CIC> {
    public:
      explicit
      PoissonChi2(CIC& cic, const Attenuation& att, const double lnSCut)
        : CICFitter<CIC>(cic, att, cic.GetNBins() - 1), fLnSCut(lnSCut) { }

      virtual
      double
      operator()(const std::vector<double>& pars)
        const
        override
      {
        const auto att = [this, &pars](const double sin2theta) {
          return this->GetAttenuation(sin2theta, pars);
        };
        return this->fCIC->GetPoissonChi2(att, fLnSCut);
      }

    private:
      double fLnSCut = 0;
    };


    template<class CIC, class Binning>
    class PoissonChi22D : public CICFitter<CIC> {
    public:
      explicit
      PoissonChi22D(CIC& cic, const Attenuation& att, const Binning& lnSBinning) :
        CICFitter<CIC>(cic, att, (cic.GetNBins() - 1)*lnSBinning.GetNBins()),
        fLnSBinning(lnSBinning)
      { }

      virtual
      double
      operator()(const std::vector<double>& pars)
        const
        override
      {
        const auto att = [this, &pars](const double sin2theta) {
          return this->GetAttenuation(sin2theta, pars);
        };
        return this->fCIC->GetPoissonChi22D(att, fLnSBinning);
      }

    private:
      Binning fLnSBinning;
    };


    template<class CIC, class Binning>
    class AndersonDarling2D : public CICFitter<CIC> {
    public:
      explicit
      AndersonDarling2D(CIC& cic, const Attenuation& att, const Binning& lnSBinning) :
        CICFitter<CIC>(cic, att, lnSBinning.GetNBins()),
        fLnSBinning(lnSBinning)
      { }

      virtual
      double
      operator()(const std::vector<double>& pars)
        const
        override
      {
        const auto att = [this, &pars](const double sin2theta) {
          return this->GetAttenuation(sin2theta, pars);
        };
        return this->fCIC->GetAndersonDarling2D(att, fLnSBinning);
      }

    private:
      Binning fLnSBinning;
    };

  }

}

#endif