TemplateFitter.cxx

Go to the documentation of this file.

//***************************************************************************************
//***************************************************************************************
 
#include "NAGASH/TemplateFitter.h"
 
using namespace NAGASH;
 
TemplateFitter::TemplateFitterFCN::TemplateFitterFCN(TH1D *_htarget, const std::vector<TH1D *> &vh, int _min_, int _max_)
{
    h_target = _htarget;
    h_template = vh;
    min = _min_;
    max = _max_;
}
 
void TemplateFitter::TemplateFitterFCN::SetRelation(std::function<double(const std::vector<double> &, const std::vector<double> &)> cf,
                                                    std::function<double(const std::vector<double> &, const std::vector<double> &, const std::vector<double> &)> ef)
{
    contentfunction = cf;
    errorfunction = ef;
}
 
double TemplateFitter::TemplateFitterFCN::operator()(const std::vector<double> &vpara) const
{
    double chi2 = 0;
 
    std::vector<double> contents(h_template.size());
    std::vector<double> errors(h_template.size());
 
    bool use_user_defined_function = bool(contentfunction) && bool(contentfunction);
 
    for (int i = min; i <= max; i++)
    {
        double BinContent = h_target->GetBinContent(i);
        double BinError = pow(h_target->GetBinError(i), 2);
 
        for (size_t j = 0; j < h_template.size(); j++)
        {
            contents[j] = h_template[j]->GetBinContent(i);
            errors[j] = h_template[j]->GetBinError(i);
        }
 
        if (use_user_defined_function)
        {
            BinContent -= contentfunction(vpara, contents);
            BinError += errorfunction(vpara, contents, errors);
        }
        else
        {
            for (size_t j = 0; j < h_template.size(); j++)
            {
                BinContent -= vpara[j] * h_template[j]->GetBinContent(i);
                BinError += pow(vpara[j] * h_template[j]->GetBinError(i), 2);
            }
        }
 
        if (BinError != 0) // avoid nan passed to chi2
            chi2 += BinContent * BinContent / BinError;
    }
 
    return chi2;
}
 
void TemplateFitter::Clear()
{
    ht = nullptr;
    vh.clear();
    sf.clear();
    sferr.clear();
    isfixed.clear();
    chi2 = 0;
}
 
void TemplateFitter::SetTarget(TH1D *h)
{
    if (!h)
    {
        MSGUser()->StartTitle("TemplateFitter::SetTarget");
        MSGUser()->MSG(MSGLevel::ERROR, "Target hist is NULL");
        MSGUser()->EndTitle();
        return;
    }
 
    ht = h;
    FitRangeMin = 1;
    FitRangeMax = ht->GetNbinsX();
}
 
void TemplateFitter::SetTemplate(TH1D *h, double SF, bool fix)
{
    if (!h)
    {
        MSGUser()->StartTitle("TemplateFitter::SetTemplate");
        MSGUser()->MSG(MSGLevel::ERROR, "Template hist is NULL");
        MSGUser()->EndTitle();
        return;
    }
 
    vh.push_back(h);
    sf.push_back(SF);
    sfmin.push_back(-999);
    sfmax.push_back(-1000);
    sferr.push_back(0);
    isfixed.push_back(fix);
}
 
void TemplateFitter::SetTemplate(TH1D *h, double SF, double SFMIN, double SFMAX)
{
    auto titleguard = MSGUser()->StartTitleWithGuard("TemplateFitter::SetTemplate");
    if (!h)
    {
        MSGUser()->MSG(MSGLevel::ERROR, "Template hist is NULL");
        return;
    }
 
    vh.push_back(h);
    sf.push_back(SF);
    if (SFMIN < SFMAX && SF < SFMAX && SF > SFMIN)
    {
        sfmin.push_back(SFMIN);
        sfmax.push_back(SFMAX);
    }
    else
    {
        MSGUser()->MSG_WARNING("Input scale factor range is wrong, brought back to default");
        sfmin.push_back(-999);
        sfmax.push_back(-1000);
    }
    sferr.push_back(0);
    isfixed.push_back(false);
}
 
void TemplateFitter::SetRelation(std::function<double(const std::vector<double> &, const std::vector<double> &)> cf,
                                 std::function<double(const std::vector<double> &, const std::vector<double> &, const std::vector<double> &)> ef)
{
    this->contentfunction = cf;
    this->errorfunction = ef;
}
 
void TemplateFitter::SetFitRange(int min, int max)
{
    if (min <= max && min > 0)
    {
        FitRangeMin = min;
        FitRangeMax = max;
    }
    else
    {
        MSGUser()->StartTitle("TemplateFitter::SetFitRange");
        MSGUser()->MSG(MSGLevel::ERROR, "Wrong fit range, brought back to default");
        MSGUser()->EndTitle();
    }
}
 
bool TemplateFitter::Check()
{
    auto titleguard = MSGUser()->StartTitleWithGuard("TemplateFitter::Check");
    if (!ht)
    {
        MSGUser()->MSG(MSGLevel::ERROR, "Target hist not set");
        return false;
    }
 
    if (sf.size() < 1)
    {
        MSGUser()->MSG(MSGLevel::ERROR, "No template hist");
        return false;
    }
 
    for (size_t i = 0; i < vh.size(); i++)
    {
        if (vh[i]->GetNbinsX() < FitRangeMax)
        {
            MSGUser()->MSG(MSGLevel::ERROR, "Wrong fit range or wrong input template hist");
            return false;
        }
    }
 
    return true;
}
 
void TemplateFitter::Fit()
{
    auto titleguard = MSGUser()->StartTitleWithGuard("TemplateFitter::Fit");
    if (!Check())
    {
        MSGUser()->MSG(MSGLevel::WARNING, "Fit won't proceed");
        return;
    }
 
    TemplateFitterFCN fcn(ht, vh, FitRangeMin, FitRangeMax);
    if (contentfunction && errorfunction)
        fcn.SetRelation(contentfunction, errorfunction);
 
    ROOT::Minuit2::MnUserParameters upar;
    for (size_t i = 0; i < sf.size(); i++)
    {
        TString name;
        name.Form("sf_%lu", i);
        if (sfmin[i] > sfmax[i])
        {
            upar.Add(name.Data(), sf[i], fabs(sf[i] / 5));
            // upar.SetLowerLimit(i, 0);
        }
        else
        {
            upar.Add(name.Data(), sf[i], (sfmax[i] - sfmin[i]) / 5, sfmin[i], sfmax[i]);
        }
        if (isfixed[i])
            upar.Fix(i);
    }
 
    // use simplex first
    ROOT::Minuit2::MnSimplex simplex(fcn, upar, 2);
    auto minsimplex = simplex(1e6);
    // then use migrad
    ROOT::Minuit2::MnMigrad migrad(fcn, minsimplex.UserState(), ROOT::Minuit2::MnStrategy(2));
    auto minmigrad = migrad(1e6);
 
    for (size_t i = 0; i < sf.size(); i++)
    {
        sf[i] = minmigrad.UserState().Value(i);
        sferr[i] = minmigrad.UserState().Error(i);
    }
    chi2 = minmigrad.Fval();
}
 
double TemplateFitter::GetSF(uint64_t i)
{
    if (i < sf.size() && i >= 0)
    {
        return sf[i];
    }
    else
    {
        MSGUser()->StartTitle("TemplateFitter::GetSF");
        MSGUser()->MSG(MSGLevel::ERROR, "No SF correspond to input number, return -1");
        MSGUser()->EndTitle();
        return -1;
    }
}
 
double TemplateFitter::GetSFError(uint64_t i)
{
    if (i < sf.size() && i >= 0)
    {
        return sferr[i];
    }
    else
    {
        MSGUser()->StartTitle("TemplateFitter::GetSFError");
        MSGUser()->MSG(MSGLevel::ERROR, "No SF correspond to input number, return -1");
        MSGUser()->EndTitle();
        return -1;
    }
}
 
double TemplateFitter::GetChi2()
{
    return chi2;
}