UnfoldTool.cxx

Go to the documentation of this file.

#include "NAGASH/UnfoldTool.h"
 
using namespace NAGASH;
 
UnfoldTool::UnfoldTool(std::shared_ptr<MSGTool> msg, int num)
    : Tool(msg)
{
    binnum = num;
 
    Data_Reco = std::vector<double>(binnum, 0);
    Data_Reco_Correct = std::vector<double>(binnum, 0);
    Theory_Truth = std::vector<double>(binnum, 0);
    Theory_Truth_Correct = std::vector<double>(binnum, 0);
    Data_Truth = std::vector<double>(binnum, 0);
    Data_Truth_Correct = std::vector<double>(binnum, 0);
    Data_Reco_Estimated = std::vector<double>(binnum, 0);
    Data_Truth_Estimated = std::vector<double>(binnum, 0);
 
    IsMeasured = std::vector<bool>(binnum, false);
    IsSet = std::vector<bool>(binnum, false);
 
    R_Matrix = std::vector<std::vector<double>>(binnum, std::vector<double>(binnum, 0));
    M_Matrix = std::vector<std::vector<double>>(binnum, std::vector<double>(binnum, 0));
    MatrixIsSet = std::vector<std::vector<bool>>(binnum, std::vector<bool>(binnum, 0));
 
    FiducialCorrection = std::vector<double>(binnum, 0);
    EfficiencyCorrection = std::vector<double>(binnum, 0);
    FCorrectionIsSet = std::vector<bool>(binnum, false);
    ECorrectionIsSet = std::vector<bool>(binnum, false);
 
    ReweightFactor = std::vector<double>(binnum, 0);
 
    Data_Truth_Correct_EachItr = std::vector<std::vector<double>>(binnum);
 
    for (int i = 0; i < binnum; i++)
        for (int j = 0; j < binnum; j++)
            MatrixIsSet[i][j] = false;
}
 
void UnfoldTool::Clear()
{
    for (int i = 0; i < binnum; i++)
    {
        IsMeasured[i] = false;
        IsSet[i] = false;
        FCorrectionIsSet[i] = false;
        ECorrectionIsSet[i] = false;
        for (int j = 0; j < binnum; j++)
            MatrixIsSet[i][j] = false;
        Data_Truth_Correct_EachItr.clear();
    }
}
 
void UnfoldTool::SetData(int index, double value)
{
    auto guard = MSGUser()->StartTitleWithGuard("UnfoldTool::SetData");
    if (index < 0)
    {
        MSGUser()->MSG_ERROR("Negative index!!! Index should be positive!!!");
        return;
    }
 
    if (index >= binnum)
    {
        MSGUser()->MSG_ERROR("Index larger than binnum");
        return;
    }
 
    MSGUser()->MSG_DEBUG("The value of No.", index, " bin for Reco-Data is set to ", value);
    Data_Reco[index] = value;
    IsMeasured[index] = true;
}
 
void UnfoldTool::SetTheory(int index, double value)
{
    auto guard = MSGUser()->StartTitleWithGuard("UnfoldTool::SetTheory");
    if (index < 0)
    {
        MSGUser()->MSG_ERROR("Negative index!!! Index should be positive!!!");
        return;
    }
 
    if (index >= binnum)
    {
        MSGUser()->MSG_ERROR("Index larger than binnum");
        return;
    }
 
    MSGUser()->MSG_DEBUG("The value of No.", index, " bin for Truth-MC is set to ", value);
    Theory_Truth[index] = value;
    IsSet[index] = true;
}
 
void UnfoldTool::SetMatrix(int indexi, int indexj, double value)
{
    auto guard = MSGUser()->StartTitleWithGuard("UnfoldTool::SetMatrix");
    if (indexi < 0)
    {
        MSGUser()->MSG_ERROR("Negative index!!! Index should be positive!!!");
        return;
    }
 
    if (indexi >= binnum)
    {
        MSGUser()->MSG_ERROR("Index larger than binnum");
        return;
    }
 
    if (indexj < 0)
    {
        MSGUser()->MSG_ERROR("Negative index!!! Index should be positive!!!");
        return;
    }
 
    if (indexj >= binnum)
    {
        MSGUser()->MSG_ERROR("Index larger than binnum");
        return;
    }
 
    MSGUser()->MSG_DEBUG("The value of No.", indexi, ",", indexj, " bin for R-Matrix is set to ", value);
    R_Matrix[indexi][indexj] = value;
    MatrixIsSet[indexi][indexj] = true;
}
 
void UnfoldTool::SetFiducialCorrection(int index, double value)
{
    auto guard = MSGUser()->StartTitleWithGuard("UnfoldTool::SetFiducial");
    if (index < 0)
    {
        MSGUser()->MSG_ERROR("Negative index!!! Index should be positive!!!");
        return;
    }
 
    if (index >= binnum)
    {
        MSGUser()->MSG_ERROR("Index larger than binnum");
        return;
    }
 
    MSGUser()->MSG_DEBUG("The value of No.", index, " bin for Fiducial is set to ", value);
    FiducialCorrection[index] = value;
    FCorrectionIsSet[index] = true;
}
 
void UnfoldTool::SetEfficiencyCorrection(int index, double value)
{
    auto guard = MSGUser()->StartTitleWithGuard("UnfoldTool::SetEfficiency");
    if (index < 0)
    {
        MSGUser()->MSG_ERROR("Negative index!!! Index should be positive!!!");
        return;
    }
 
    if (index >= binnum)
    {
        MSGUser()->MSG_ERROR("Index larger than binnum");
        return;
    }
 
    MSGUser()->MSG_DEBUG("The value of No.", index, " bin for Efficiency is set to ", value);
 
    EfficiencyCorrection[index] = value;
    ECorrectionIsSet[index] = true;
}
 
void UnfoldTool::BayesMatrix()
{
    // Mij = Rji * Truthi / Recoj = Rji * Truthi / Sigmak(Rjk* Truthk)
    for (int j = 0; j < binnum; j++)
    {
        double down = 0;
        for (int k = 0; k < binnum; k++)
            down = down + R_Matrix[j][k] * Theory_Truth_Correct[k];
 
        for (int i = 0; i < binnum; i++)
            M_Matrix[i][j] = R_Matrix[j][i] * Theory_Truth_Correct[i] / down;
    }
}
 
void UnfoldTool::EstimateDataReco()
{
    for (int i = 0; i < binnum; i++)
    {
        if (IsMeasured[i])
            continue;
 
        Data_Reco_Estimated[i] = 0;
        for (int j = 0; j < binnum; j++)
            Data_Reco_Estimated[i] = Data_Reco_Estimated[i] + R_Matrix[i][j] * Theory_Truth_Correct[j];
    }
}
 
void UnfoldTool::EstimateDataTruth()
{
    for (int i = 0; i < binnum; i++)
    {
        if (IsMeasured[i])
        {
            Data_Truth_Estimated[i] = 0;
            for (int j = 0; j < binnum; j++)
                if (IsMeasured[j])
                    Data_Truth_Estimated[i] = Data_Truth_Estimated[i] + M_Matrix[i][j] * Data_Reco_Correct[j];
                else
                    Data_Truth_Estimated[i] = Data_Truth_Estimated[i] + M_Matrix[i][j] * Data_Reco_Estimated[j];
 
            Data_Truth[i] = Data_Truth_Estimated[i];
            Data_Truth_Correct[i] = Data_Truth[i] / EfficiencyCorrection[i];
        }
    }
}
 
void UnfoldTool::RenewTheory()
{
    for (int i = 0; i < binnum; i++)
    {
        if (IsMeasured[i])
            Theory_Truth_Correct[i] = Data_Truth_Estimated[i];
    }
}
 
void UnfoldTool::PrintResult()
{
    auto guard = MSGUser()->StartTitleWithGuard("UnfoldTool::PrintResult");
    MSGUser()->MSG_DEBUG("Binnum   Data   Theory   Original-Data");
    for (int i = 0; i < binnum; i++)
        if (IsMeasured[i])
            MSGUser()->MSG_DEBUG(i, "      ", Data_Truth_Correct[i], "   ", Theory_Truth[i], "   ", Data_Reco[i]);
        else
            MSGUser()->MSG_DEBUG(i, "      Fixed to Theory   ", Theory_Truth[i], "   ", Data_Reco_Estimated[i] / FiducialCorrection[i]);
 
    for (int i = 0; i < binnum; i++)
        if (IsMeasured[i])
            Data_Truth_Correct_EachItr[i].emplace_back(Data_Truth_Correct[i]);
        else
            Data_Truth_Correct_EachItr[i].emplace_back(Theory_Truth[i]);
}
 
void UnfoldTool::BayesUnfold(int iteration)
{
    auto guard = MSGUser()->StartTitleWithGuard("UnfoldTool::BayesUnfold");
    if (!Check())
    {
        MSGUser()->MSG_ERROR("Check Failed");
        return;
    }
 
    for (int i = 0; i < binnum; i++)
    {
        if (IsMeasured[i])
            Data_Reco_Correct[i] = Data_Reco[i] * FiducialCorrection[i];
        Theory_Truth_Correct[i] = Theory_Truth[i] * EfficiencyCorrection[i];
    }
 
    bool SaveLess = false;
    if (iteration > 10000)
    {
        SaveLess = true;
        MSGUser()->MSG_DEBUG("Too many iterations, save only 1, 10, 100, ... iteration");
    }
 
    int saveitr = 1;
    for (int i = 0; i < iteration; i++)
    {
        MSGUser()->MSG_DEBUG("Iteration ", i, " Start");
        BayesMatrix();
        EstimateDataReco();
        EstimateDataTruth();
        RenewTheory();
 
        if (i == saveitr && SaveLess)
        {
            PrintResult();
            saveitr = saveitr * 10;
        }
        else if (!SaveLess)
            PrintResult();
 
        MSGUser()->MSG_DEBUG("Iteration ", i, " End");
    }
}
 
void UnfoldTool::EstimateReweightFactor()
{
    for (int i = 0; i < binnum; i++)
    {
        Data_Reco_Estimated[i] = 0;
        for (int j = 0; j < binnum; j++)
            Data_Reco_Estimated[i] = Data_Reco_Estimated[i] + R_Matrix[i][j] * Theory_Truth_Correct[j];
    }
 
    for (int i = 0; i < binnum; i++)
        ReweightFactor[i] = Data_Reco_Correct[i] / Data_Reco_Estimated[i];
}
 
void UnfoldTool::ApplyReweightFactor()
{
    for (int i = 0; i < binnum; i++)
    {
        Data_Truth[i] = ReweightFactor[i] * Theory_Truth_Correct[i];
        Data_Truth_Estimated[i] = ReweightFactor[i] * Theory_Truth_Correct[i];
        Data_Truth_Correct[i] = Data_Truth[i] / EfficiencyCorrection[i];
    }
}
 
void UnfoldTool::ReweightUnfold(int &iteration)
{
    auto guard = MSGUser()->StartTitleWithGuard("UnfoldTool::ReweightUnfold");
    if (!Check())
    {
        MSGUser()->MSG_ERROR("Check Failed");
        return;
    }
 
    // Initialize ReweightFactor
    for (int i = 0; i < binnum; i++)
        ReweightFactor[i] = 1;
 
    for (int i = 0; i < binnum; i++)
    {
        if (IsMeasured[i])
            Data_Reco_Correct[i] = Data_Reco[i] * FiducialCorrection[i];
        Theory_Truth_Correct[i] = Theory_Truth[i] * EfficiencyCorrection[i];
    }
 
    bool keepdoing = false;
    int reweight_count = 0;
    double chi2 = 0;
    do
    {
        keepdoing = false;
        EstimateReweightFactor();
        ApplyReweightFactor();
 
        reweight_count++;
 
        chi2 = 0;
        for (int j = 0; j < binnum; j++)
        {
            chi2 = chi2 + (Data_Reco[j] - Data_Reco_Estimated[j] / FiducialCorrection[j]) * (Data_Reco[j] - Data_Reco_Estimated[j] / FiducialCorrection[j]) / Data_Reco[j];
            if ((Data_Reco[j] - Data_Reco_Estimated[j] / FiducialCorrection[j]) * (Data_Reco[j] - Data_Reco_Estimated[j] / FiducialCorrection[j]) / Data_Reco[j] > 0.1)
                keepdoing = true;
            MSGUser()->MSG_DEBUG(Data_Reco[j], "   ",
                                 Data_Reco_Estimated[j] / FiducialCorrection[j], "   ",
                                 Data_Reco_Correct[j], "   ", Data_Reco_Estimated[j], "   ",
                                 ReweightFactor[j], "   ",
                                 Theory_Truth_Correct[j], "   ",
                                 Theory_Truth_Correct[j] * ReweightFactor[j], "   ",
                                 (Data_Reco[j] - Data_Reco_Estimated[j] / FiducialCorrection[j]) * (Data_Reco[j] - Data_Reco_Estimated[j] / FiducialCorrection[j]) / Data_Reco[j]);
 
            if (reweight_count > 2000)
                keepdoing = false;
        }
        MSGUser()->MSG_DEBUG("Reweight No.", reweight_count, "   Chi-square: ", chi2);
        RenewTheory();
    } while (keepdoing);
 
    MSGUser()->MSG_DEBUG("Reweight No.", reweight_count, "   Chi-square: ", chi2);
    PrintResult();
    iteration = reweight_count;
}
 
bool UnfoldTool::Check()
{
    // Check every thing need is set
    for (int i = 0; i < binnum; i++)
        if (IsSet[i] == false || FCorrectionIsSet[i] == false || ECorrectionIsSet[i] == false)
            return false;
 
    for (int i = 0; i < binnum; i++)
        for (int j = 0; j < binnum; j++)
            if (MatrixIsSet[i][j] == false)
                return false;
 
    return true;
}
 
void UnfoldTool::DisplayEachItr()
{
    for (int i = 0; i < Data_Truth_Correct_EachItr[0].size(); i++)
    {
        std::cout << "Iteration No." << i << "    ";
        for (int j = 0; j < binnum; j++)
            std::cout << Data_Truth_Correct_EachItr[j][i] << "  ";
        std::cout << std::endl;
    }
}
 
void UnfoldTool::GetResult(int itration, TH1D *result)
{
    auto guard = MSGUser()->StartTitleWithGuard("UnfoldTool::GetResult");
    if (result == nullptr)
    {
        MSGUser()->MSG_ERROR("Hist for result is not set!!!");
        return;
    }
 
    if (result->GetNbinsX() < binnum)
    {
        MSGUser()->MSG_ERROR("Not enought bins for the result!!!");
        return;
    }
 
    for (int i = 0; i < binnum; i++)
    {
        result->SetBinContent(i + 1, Data_Truth_Correct_EachItr[i][itration]);
        result->SetBinError(i + 1, 0);
    }
}