Chi2Fitter.cxx

Go to the documentation of this file.

//**************************************************************************************
//**************************************************************************************
 
#include "NAGASH/Chi2Fitter.h"
 
using namespace NAGASH;
using namespace std;
 
Chi2Fitter::Chi2Fitter(std::shared_ptr<MSGTool> msg, uint64_t dim, TH1D *_target_, TH2D *_cov_)
    : Tool(msg), dimension(dim), target(_target_), targetcov(_cov_), cov(dim, dim), histtool(msg)
{
    auto titleguard = MSGUser()->StartTitleWithGuard("Chi2Fitter::Chi2Fitter");
    if (!target)
    {
        MSGUser()->MSG_ERROR("Target histogram is not valid");
    }
    else
    {
        RangeMin = 1;
        RangeMax = target->GetNbinsX();
        isvalid = true;
    }
}
 
Chi2Fitter::Chi2Fitter(std::shared_ptr<MSGTool> msg, uint64_t dim)
    : Tool(msg), dimension(dim), cov(dim, dim), histtool(msg)
{
    isvalid = true;
}
 
void Chi2Fitter::SetRange(int min, int max)
{
    auto titleguard = MSGUser()->StartTitleWithGuard("Chi2Fitter::SetRange");
    if (min > max || min <= 0 || max > RangeMax)
    {
        MSGUser()->MSG_WARNING("Wrong range input, brought back to default");
    }
    else
    {
        RangeMin = min;
        RangeMax = max;
    }
}
 
void Chi2Fitter::SetModel(const std::vector<double> &vv, double _chi2_)
{
    auto titleguard = MSGUser()->StartTitleWithGuard("Chi2Fitter::SetModel");
    if (vv.size() != dimension)
    {
        MSGUser()->MSG_ERROR("Size of list of model parameters does not match");
        return;
    }
 
    vchi2.emplace_back(_chi2_);
    vvalue.emplace_back(vv);
 
    // prepare for normalization
    if (vchi2.size() == 1)
    {
        vvaluemax = vv;
        vvaluemin = vv;
    }
    else
    {
        for (size_t i = 0; i < dimension; i++)
        {
            if (vv[i] > vvaluemax[i])
                vvaluemax[i] = vv[i];
            if (vv[i] < vvaluemin[i])
                vvaluemin[i] = vv[i];
        }
    }
}
 
void Chi2Fitter::SetModel(const std::vector<double> &vv, TH1D *hmodel, TH2D *covmodel)
{
    if (!isvalid)
        return;
 
    auto titleguard = MSGUser()->StartTitleWithGuard("Chi2Fitter::SetModel");
    if (!hmodel)
    {
        MSGUser()->MSG_ERROR("Model histogram is not valid");
        return;
    }
    if (!target)
    {
        MSGUser()->MSG_ERROR("Target histogram is not valid");
        return;
    }
 
    if (hmodel->GetNbinsX() < RangeMax)
    {
        MSGUser()->MSG_ERROR("Model histogram has less bins than fit range");
        return;
    }
 
    if (covmodel && (hmodel->GetNbinsX() != covmodel->GetNbinsX() || hmodel->GetNbinsX() != covmodel->GetNbinsY()))
    {
        MSGUser()->MSG_ERROR("Model histogram and covariance matrix do not match");
        return;
    }
 
    if (vv.size() != dimension)
    {
        MSGUser()->MSG_ERROR("Size of list of model parameters does not match");
        return;
    }
 
    double chi2 = 0;
    if (!covmodel && !targetcov)
    {
        // no cov
        for (int i = RangeMin; i <= RangeMax; i++)
        {
            if (target->GetBinError(i) != 0 || hmodel->GetBinError(i) != 0)
                chi2 += pow(target->GetBinContent(i) - hmodel->GetBinContent(i), 2) / (target->GetBinError(i) * target->GetBinError(i) + hmodel->GetBinError(i) * hmodel->GetBinError(i));
        }
    }
    else
    {
        TMatrixD matcovtarget(target->GetNbinsX(), target->GetNbinsX());
        TMatrixD matcovmodel(hmodel->GetNbinsX(), hmodel->GetNbinsX());
        matcovtarget *= 0;
        matcovmodel *= 0;
 
        if (targetcov)
            matcovtarget = histtool.ConvertTH2DToMatrix(targetcov);
        else
            for (int i = 1; i <= target->GetNbinsX(); i++)
                matcovtarget(i - 1, i - 1) = target->GetBinError(i) * target->GetBinError(i);
 
        if (covmodel)
            matcovmodel = histtool.ConvertTH2DToMatrix(covmodel);
        else
            for (int i = 1; i <= hmodel->GetNbinsX(); i++)
                matcovmodel(i - 1, i - 1) = hmodel->GetBinError(i) * hmodel->GetBinError(i);
 
        int psize = RangeMax - RangeMin + 1;
        TMatrixD partialmatcov(psize, psize);
 
        for (int i = 0; i < psize; i++)
            for (int j = 0; j < psize; j++)
                partialmatcov(i, j) = matcovtarget(i + RangeMin - 1, j + RangeMin - 1) + matcovmodel(i + RangeMin - 1, j + RangeMin - 1);
 
        auto invertcov = partialmatcov.InvertFast();
        for (int i = 0; i < psize; i++)
            for (int j = 0; j < psize; j++)
                chi2 += invertcov(i, j) * (target->GetBinContent(i + RangeMin) - hmodel->GetBinContent(i + RangeMin)) *
                        (target->GetBinContent(j + RangeMin) - hmodel->GetBinContent(j + RangeMin));
    }
 
    vchi2.emplace_back(chi2);
    vvalue.emplace_back(vv);
 
    // prepare for normalization
    if (vchi2.size() == 1)
    {
        vvaluemax = vv;
        vvaluemin = vv;
    }
    else
    {
        for (size_t i = 0; i < dimension; i++)
        {
            if (vv[i] > vvaluemax[i])
                vvaluemax[i] = vv[i];
            if (vv[i] < vvaluemin[i])
                vvaluemin[i] = vv[i];
        }
    }
}
 
double Chi2Fitter::GetChi2(uint64_t index)
{
    if (index < vchi2.size() && index >= 0)
        return vchi2[index];
    else
        return -1;
}
 
double Chi2Fitter::FitChi2()
{
    return chi2;
}
 
double Chi2Fitter::Value(uint64_t index)
{
    if (index < value.size() && index >= 0)
        return value[index];
    else
        return -1;
}
 
double Chi2Fitter::Error(uint64_t index)
{
    if ((int)index < cov.GetNcols() && index >= 0)
        return sqrt(cov(index, index));
    else
        return -1;
}
 
double Chi2Fitter::OffSet()
{
    return this->offset;
}
 
// structure of parameters
// first (dim + 1) * dim / 2 of elements are lower triangular of inverse of cov
// e.g. matrix
//     4 1 0
//     1 2 5
//     0 5 9
// are stored as {4, 1, 2, 0, 5, 9}
// then "dim" of elements are parameters of the linear terms
// last is offset
 
TMatrixD RestoreMatrix(int dim, const Eigen::VectorXf &vpara)
{
    TMatrixD mat(dim, dim);
    int count = 0;
    for (int i = 0; i < dim; i++)
        for (int j = 0; j <= i; j++)
        {
            mat(i, j) = vpara(count);
            ++count;
        }
 
    // since the covariance matrix is symmetric
    for (int i = 0; i < dim; i++)
    {
        for (int j = i + 1; j < dim; j++)
        {
            mat(i, j) = mat(j, i);
        }
    }
 
    return mat;
}
 
TMatrixD RestoreMatrix(int dim, const std::vector<double> &vpara)
{
    TMatrixD mat(dim, dim);
    int count = 0;
    for (int i = 0; i < dim; i++)
        for (int j = 0; j <= i; j++)
        {
            mat(i, j) = vpara[count];
            ++count;
        }
 
    // since the covariance matrix is symmetric
    for (int i = 0; i < dim; i++)
    {
        for (int j = i + 1; j < dim; j++)
        {
            mat(i, j) = mat(j, i);
        }
    }
 
    return mat;
}
 
void Chi2Fitter::Fit(FitMethod fm)
{
    if (!isvalid)
        return;
 
    auto titleguard = MSGUser()->StartTitleWithGuard("Chi2Fitter::Fit");
    // firstly, transform the range of value into [1, 3]
    // to avoid large difference between x and x^2
    std::vector<double> vnf1, vnf2;
    std::vector<std::vector<double>> vvaluescaled;
    for (size_t i = 0; i < dimension; i++)
    {
        vnf1.emplace_back((vvaluemax[i] - vvaluemin[i]) / 2);
        vnf2.emplace_back(-vvaluemax[i] / 2 + 3 * vvaluemin[i] / 2);
        if (fabs(vnf1.back()) < 1e-8)
            vnf1.back() = 1;
 
        /*
        vnf1.emplace_back(1);
        vnf2.emplace_back(0);
        */
    }
 
    std::vector<double> vpara;
    if (fm == FitMethod::Eigen)
    {
        Eigen::MatrixXf X = Eigen::MatrixXf(vchi2.size(), (dimension * dimension + dimension) / 2 + dimension + 1);
        Eigen::VectorXf Y = Eigen::VectorXf(vchi2.size());
        for (size_t n = 0; n < vchi2.size(); n++)
        {
            int count = 0;
            for (size_t i = 0; i < dimension; i++)
                for (size_t j = 0; j <= i; j++)
                {
                    X(n, count) = ((vvalue[n][i] - vnf2[i]) / vnf1[i]) * (vvalue[n][j] - vnf2[j]) / vnf1[j];
                    ++count;
                }
 
            for (size_t i = 0; i < dimension; i++)
            {
                X(n, count) = ((vvalue[n][i] - vnf2[i]) / vnf1[i]);
                ++count;
            }
 
            X(n, count) = 1;
            Y(n) = vchi2[n];
        }
 
        /*
        for (int n = 0; n < vchi2.size(); n++)
        {
            for (int i = 0; i < X.cols(); i++)
                std::cout << X(n, i) << " ";
            std::cout << Y(n) << std::endl;
        }
        */
 
        Eigen::VectorXf B = X.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(Y); // SVD
        // Eigen::VectorXf B = X.colPivHouseholderQr().solve(Y); // QR decomposition
        // std::cout << B << std::endl;
        // std::cout << (X.transpose() * X).ldlt().solve(X.transpose() * Y) << std::endl; // direct solve
 
        // convert fitted parameters to target parameters
        for (int i = 0; i < B.size(); i++)
            vpara.emplace_back(B(i));
    }
 
    if (fm == FitMethod::Minuit)
    {
        std::vector<std::vector<double>> scaledvvalue;
        for (size_t i = 0; i < vchi2.size(); i++)
        {
            auto sv = vvalue[i];
            for (size_t j = 0; j < dimension; j++)
                sv[j] = (sv[j] - vnf2[j]) / vnf1[j];
 
            scaledvvalue.emplace_back(sv);
        }
        Chi2FitterFCN fcn(dimension, &scaledvvalue, &vchi2);
 
        ROOT::Minuit2::MnUserParameters upar;
        for (size_t i = 0; i < dimension; i++)
            for (size_t j = 0; j <= i; j++)
                upar.Add(TString::Format("invertcov_%lu_%lu", i, j).Data(), 0, 1);
 
        for (size_t i = 0; i < dimension; i++)
            upar.Add(TString::Format("linearpar_%lu", i).Data(), 2, 2);
 
        upar.Add("offset", 0, 1);
 
        // use simplex first
        ROOT::Minuit2::MnSimplex simplex(fcn, upar, 2);
        auto minsimplex = simplex(1e5);
        // then use migrad
        ROOT::Minuit2::MnMigrad migrad(fcn, minsimplex.UserState(), ROOT::Minuit2::MnStrategy(2));
        auto minmigrad = migrad(1e5);
 
        chi2 = minmigrad.Fval();
        vpara = minmigrad.UserState().Params();
    }
 
    cov = RestoreMatrix(dimension, vpara); // actually now the cov is the inverse of the covariance
    for (size_t i = 0; i < dimension; i++)
        for (size_t j = 0; j < dimension; j++)
            cov(i, j) /= vnf1[i] * vnf1[j];
 
    auto invertcov = cov;
    cov.Invert();
 
    std::vector<double> linearterm;
    value = std::vector<double>(dimension, 0);
    for (size_t i = 0; i < dimension; i++)
        linearterm.push_back(vpara[i + (dimension * dimension + dimension) / 2]);
 
    for (size_t i = 0; i < dimension; i++)
    {
        for (size_t j = 0; j < dimension; j++)
            value[i] -= cov(i, j) * linearterm[j] / (2 * vnf1[j]);
 
        value[i] += vnf2[i];
    }
 
    // calculation of the offset
    offset = vpara.back();
    auto tempcov = RestoreMatrix(dimension, vpara).Invert();
    auto tempinvertcov = RestoreMatrix(dimension, vpara);
    std::vector<double> tempvalue(dimension, 0);
    for (size_t i = 0; i < dimension; i++)
        for (size_t j = 0; j < dimension; j++)
            tempvalue[i] -= tempcov(i, j) * linearterm[j] / 2;
 
    for (size_t i = 0; i < dimension; i++)
        for (size_t j = 0; j < dimension; j++)
            offset -= tempvalue[i] * tempinvertcov(i, j) * tempvalue[j];
 
    // check the variance
    for (size_t i = 0; i < dimension; i++)
        if (cov(i, i) <= 0 || !isfinite(cov(i, i)))
            MSGUser()->MSG_WARNING("Variance of value ", i, " is negative or is not finite, the fit results are not trustworthy");
}
 
double Chi2Fitter::Chi2FitterFCN::QuadraticChi2(const std::vector<double> &vx, const std::vector<double> &vpara) const
{
    double chi2 = vpara.back();
    auto mat = RestoreMatrix(dim, vpara);
 
    for (size_t i = 0; i < dim; i++)
        for (size_t j = 0; j < dim; j++)
            chi2 += vx[i] * mat(i, j) * vx[j];
 
    for (size_t i = 0; i < dim; i++)
        chi2 += vx[i] * vpara[i + (dim * dim + dim) / 2];
 
    return chi2;
}
 
std::vector<double> Chi2Fitter::Chi2FitterFCN::GradientQuadraticChi2(const std::vector<double> &vx, const std::vector<double> &vpara) const
{
    std::vector<double> gradient(vpara.size(), 0);
    gradient.back() = 1;
    auto mat = RestoreMatrix(dim, vpara);
 
    int count = 0;
    for (size_t i = 0; i < dim; i++)
        for (size_t j = 0; j <= i; j++)
        {
            if (i == j)
                gradient[count] = vx[i] * vx[i];
            else
                gradient[count] = 2 * vx[i] * vx[j];
 
            ++count;
        }
 
    for (size_t i = 0; i < dim; i++)
        gradient[count + i] = vx[i];
 
    return gradient;
}
 
double Chi2Fitter::Chi2FitterFCN::operator()(const std::vector<double> &vpara) const
{
    double chi2 = 0;
    for (size_t i = 0; i < vc->size(); i++)
        chi2 += pow(QuadraticChi2(vv->at(i), vpara) - vc->at(i), 2);
 
    return chi2;
}
 
std::vector<double> Chi2Fitter::Chi2FitterFCN::Gradient(const std::vector<double> &vpara) const
{
    std::vector<double> gradient(vpara.size(), 0);
    for (size_t i = 0; i < vc->size(); i++)
    {
        auto chi2gradient = GradientQuadraticChi2(vv->at(i), vpara);
        double dchi2 = QuadraticChi2(vv->at(i), vpara) - vc->at(i);
        for (size_t j = 0; j < gradient.size(); j++)
            gradient[j] += 2 * dchi2 * chi2gradient[j];
    }
 
    return gradient;
}