ProfileFitter.cxx

Go to the documentation of this file.

//***************************************************************************************
//***************************************************************************************
 
#include "NAGASH/ProfileFitter.h"
#include "NAGASH/ThreadPool.h"
 
using namespace NAGASH;
using namespace std;
 
ProfileFitter::ProfileFitter(std::shared_ptr<MSGTool> msg, int nthread) : Tool(msg)
{
    NTHREAD = nthread;
    if (NTHREAD > 1)
        pool = std::make_shared<ThreadPool>(NTHREAD);
 
    vector<TString> Sample_Name_vec;
    Gamma_Map.emplace(std::pair<TString, vector<TString>>("Gamma", Sample_Name_vec));
}
 
void ProfileFitter::BookParameter(TString name, double init, double min, double max, ParameterType type)
{
    auto msgguard = MSGUser()->StartTitleWithGuard("ProfileFitter::BookParameter");
 
    Parameter par;
    par.Name = name;
    par.Init = init;
    par.Value = init;
    par.Max = max;
    par.Min = min;
    par.Type = type;
 
    if (par.Type == ParameterType::Interest)
        par.SetStrategy(POI_DEFAULT_STRATEGY);
    else
        par.SetStrategy(NUISANCE_DEFAULT_STRATEGY);
 
    auto findpar = Par_Map.find(name);
    if (findpar != Par_Map.end())
    {
        MSGUser()->MSG_ERROR("Parameter ", name.TString::Data(), " has already been defined !");
        MSGUser()->MSG_ERROR("This defination is ignored !");
        return;
    }
    else
    {
        Par_Map.emplace(std::pair<TString, Parameter>(name, par));
    }
}
 
void ProfileFitter::BookObservableData(TString name, const std::vector<double> &vec, const std::vector<double> &vec_stat)
{
    auto msgguard = MSGUser()->StartTitleWithGuard("ProfileFitter::BookObservableData");
 
    auto findobs = Obs_Map.find(name);
    Observable *obs = nullptr;
    if (findobs != Obs_Map.end())
    {
        obs = &(findobs->second);
    }
    else
    {
        obs = InitializeNewObservable(name);
    }
 
    if (obs->isdata_booked)
    {
        MSGUser()->MSG_ERROR("The data of Observable ", name.TString::Data(), " has already been booked !");
        MSGUser()->MSG_ERROR("This new booking is ignored !");
        return;
    }
 
    obs->Data = vec;
    obs->Profile = vec;
    obs->Data_Stat = vec_stat;
    obs->isdata_booked = true;
}
 
void ProfileFitter::BookObservableNominal(TString name, TString samplename, const std::vector<double> &vec, const std::vector<double> &vec_stat)
{
    auto msgguard = MSGUser()->StartTitleWithGuard("ProfileFitter::BookObservableNominal");
 
    bool isfindgamma = false;
    for (auto &g : Gamma_Map)
    {
        if (g.second.size() > 0)
        {
            for (int i = 0; i < g.second.size(); i++)
            {
                if (g.second[i].TString::EqualTo(samplename.TString::Data()))
                    isfindgamma = true;
            }
        }
    }
    if (!isfindgamma)
        Gamma_Map["Gamma"].emplace_back(samplename);
 
    auto findobs = Obs_Map.find(name);
    Observable *obs = nullptr;
    if (findobs != Obs_Map.end())
    {
        obs = &(findobs->second);
    }
    else
    {
        obs = InitializeNewObservable(name);
    }
 
    auto findsample = obs->Sample_Map.find(samplename);
    Observable::Sample *samp = nullptr;
    if (findsample != obs->Sample_Map.end())
    {
        samp = &(findsample->second);
        if (samp->isnominal_booked)
        {
            MSGUser()->MSG_ERROR("The nominal of Sample ", samplename.TString::Data(), " in Observable ", name.TString::Data(), " has already been booked !");
            MSGUser()->MSG_ERROR("This new booking is ignored !");
            return;
        }
    }
    else
    {
        Observable::Sample newsamp;
        newsamp.Name = samplename;
        obs->Sample_Map.emplace(std::pair<TString, ProfileFitter::Observable::Sample>(samplename, newsamp));
        samp = &((obs->Sample_Map)[samplename]);
 
        obs->N_Samples = obs->N_Samples + 1;
    }
    samp->Nominal = vec;
    samp->Profile = vec;
    samp->Nominal_Stat = vec_stat;
    samp->isnominal_booked = true;
 
    // check if any value below zero
    for (int i = 0; i < vec.size(); i++)
    {
        if (vec[i] < 0)
        {
            MSGUser()->MSG_WARNING("Nominal value of sample ", samplename, " at bin ", i, " is below zero, set to 1e-5! Orignal value and error ", vec[i], " ", vec_stat[i]);
            samp->Nominal[i] = 1e-5;
            samp->Profile[i] = 1e-5;
        }
    }
}
 
ProfileFitter::Parameter *ProfileFitter::GetParameter(TString name)
{
    auto msgguard = MSGUser()->StartTitleWithGuard("ProfileFitter::GetParamter");
 
    auto findpar = Par_Map.find(name);
    if (findpar != Par_Map.end())
        return &(findpar->second);
 
    MSGUser()->MSG_ERROR("Parameter ", name.TString::Data(), " has not been defined !");
    return nullptr;
}
 
void ProfileFitter::BookObservableVariation(TString name, TString samplename, const std::vector<double> &vec, TString parname, double sigma)
{
    auto msgguard = MSGUser()->StartTitleWithGuard("ProfileFitter::BookObservableVariation");
 
    bool isfindgamma = false;
    for (auto &g : Gamma_Map)
    {
        if (g.second.size() > 0)
        {
            for (int i = 0; i < g.second.size(); i++)
            {
                if (g.second[i].TString::EqualTo(samplename.TString::Data()))
                    isfindgamma = true;
            }
        }
    }
    if (!isfindgamma)
        Gamma_Map["Gamma"].emplace_back(samplename);
 
    auto findpar = Par_Map.find(parname);
    if (findpar == Par_Map.end())
    {
        MSGUser()->MSG_ERROR("Parameter ", parname.TString::Data(), " has not been defined !");
        MSGUser()->MSG_ERROR("This Variation has been ignored !");
        return;
    }
 
    // check the input vector
    std::vector<double> modified_vec = vec;
    for (int i = 0; i < vec.size(); i++)
    {
        if (vec[i] < 0)
        {
            modified_vec[i] = 1e-5;
            MSGUser()->MSG_WARNING("Variation value of sample ", samplename, " at bin ", i, " of par ", parname, " is below zero, set to 1e-5! Orignal value ", vec[i]);
        }
    }
 
    Parameter *par = GetParameter(parname);
    Observable *obs = nullptr;
 
    auto findobs = Obs_Map.find(name);
    if (findobs != Obs_Map.end())
    {
        obs = &(findobs->second);
    }
    else
    {
        obs = InitializeNewObservable(name);
    }
 
    auto findsample = obs->Sample_Map.find(samplename);
    Observable::Sample *samp = nullptr;
    if (findsample != obs->Sample_Map.end())
    {
        samp = &(findsample->second);
    }
    else
    {
        Observable::Sample newsamp;
        newsamp.Name = samplename;
        obs->Sample_Map.emplace(std::pair<TString, ProfileFitter::Observable::Sample>(samplename, newsamp));
        samp = &((obs->Sample_Map)[samplename]);
 
        obs->N_Samples = obs->N_Samples + 1;
    }
 
    auto findsigma = samp->Variation_Sigma_Map.find(par);
    if (findsigma != samp->Variation_Sigma_Map.end())
    {
        findsigma->second.emplace_back(sigma);
    }
    else
    {
        vector<double> sigma_vec;
        sigma_vec.emplace_back(sigma);
        samp->Variation_Sigma_Map.emplace(std::pair<Parameter *, vector<double>>(par, sigma_vec));
    }
 
    auto findvar = samp->Variation_Map.find(par);
    if (findvar != samp->Variation_Map.end())
    {
        findvar->second.emplace_back(modified_vec);
    }
    else
    {
        vector<vector<double>> var_vec;
        var_vec.emplace_back(modified_vec);
        samp->Variation_Map.emplace(std::pair<Parameter *, vector<vector<double>>>(par, var_vec));
        vector<std::shared_ptr<ProfileCache>> cache_vec;
        for (int index = 0; index < vec.size(); index++)
            cache_vec.emplace_back(nullptr);
        samp->Variation_Cache.emplace(std::pair<Parameter *, vector<std::shared_ptr<ProfileCache>>>(par, cache_vec));
    }
 
    auto findsym = samp->NeedSymmetrization_Map.find(par);
    if (findsym == samp->NeedSymmetrization_Map.end())
    {
        samp->NeedSymmetrization_Map.emplace(std::pair<Parameter *, bool>(par, false));
    }
 
    auto findsmooth = samp->SmoothLevel_Map.find(par);
    if (findsmooth == samp->SmoothLevel_Map.end())
    {
        samp->SmoothLevel_Map.emplace(std::pair<Parameter *, int>(par, -1));
    }
 
    samp->N_Variations = samp->N_Variations + 1;
 
    if (fabs(sigma) < 1e-5)
    {
        auto findvar = samp->Variation_Nominal_Map.find(par);
        if (findvar != samp->Variation_Nominal_Map.end())
        {
            findvar->second = modified_vec;
            MSGUser()->MSG_ERROR("Nominal for Parameter ", par->Name.TString::Data(), " Sample ", samplename.TString::Data(), " has already defined !");
        }
        else
        {
            samp->Variation_Nominal_Map.emplace(std::pair<Parameter *, vector<double>>(par, modified_vec));
        }
    } // 0 sigma
    if (fabs(sigma - 1) < 1e-5)
    {
        auto findvar = samp->Variation_Up_Map.find(par);
        if (findvar != samp->Variation_Up_Map.end())
        {
            findvar->second = modified_vec;
            MSGUser()->MSG_ERROR("Up Variation for Parameter ", par->Name.TString::Data(), " Sample ", samplename.TString::Data(), " has already defined !");
        }
        else
        {
            samp->Variation_Up_Map.emplace(std::pair<Parameter *, vector<double>>(par, modified_vec));
        }
    } // 1 sigma
    if (fabs(sigma + 1) < 1e-5)
    {
        auto findvar = samp->Variation_Down_Map.find(par);
        if (findvar != samp->Variation_Down_Map.end())
        {
            findvar->second = modified_vec;
            MSGUser()->MSG_ERROR("Down Variation for Parameter ", par->Name.TString::Data(), " Sample ", samplename.TString::Data(), " has already defined !");
        }
        else
        {
            samp->Variation_Down_Map.emplace(std::pair<Parameter *, vector<double>>(par, modified_vec));
        }
    } // -1 sigma
}
 
void ProfileFitter::BookObservableVariationFromUnc(TString name, TString samplename, const std::vector<double> &vec, const std::vector<double> &vec_unc, TString parname)
{
    auto msgguard = MSGUser()->StartTitleWithGuard("ProfileFitter::BookObservableVariationUnc");
 
    bool isfindgamma = false;
    for (auto &g : Gamma_Map)
    {
        if (g.second.size() > 0)
        {
            for (int i = 0; i < g.second.size(); i++)
            {
                if (g.second[i].TString::EqualTo(samplename.TString::Data()))
                    isfindgamma = true;
            }
        }
    }
    if (!isfindgamma)
        Gamma_Map["Gamma"].emplace_back(samplename);
 
    auto findpar = Par_Map.find(parname);
    if (findpar == Par_Map.end())
    {
        MSGUser()->MSG_ERROR("Parameter ", parname.TString::Data(), " has not been defined !");
        MSGUser()->MSG_ERROR("This Variation has been ignored !");
        return;
    }
 
    std::vector<double> vec_up = vec;
    std::vector<double> vec_down = vec;
 
    for (int i = 0; i < vec.size(); i++)
    {
        vec_up[i] += vec_unc[i];
        vec_down[i] -= vec_unc[i];
    }
 
    BookObservableVariationNominal(name, samplename, vec, parname);
    BookObservableVariationUp(name, samplename, vec_up, parname);
    BookObservableVariationDown(name, samplename, vec_down, parname);
}
 
ProfileFitter::Observable *ProfileFitter::InitializeNewObservable(TString name)
{
    auto msgguard = MSGUser()->StartTitleWithGuard("ProfileFitter::InitializeNewObservable");
 
    auto findobs = Obs_Map.find(name);
    if (findobs != Obs_Map.end())
    {
        MSGUser()->MSG_ERROR("Observable ", name.TString::Data(), " has already been defined !");
        return nullptr;
    }
 
    ProfileFitter::Observable obs;
    obs.Name = name;
    Obs_Map.emplace(std::pair<TString, ProfileFitter::Observable>(name, obs));
 
    return &(Obs_Map[name]);
}
 
// Main Fitting
void ProfileFitter::Fit(ProfileFitter::FitMethod method)
{
    auto msgguard = MSGUser()->StartTitleWithGuard("ProfileFitter::Fit");
    PrintInfo();
    Prepare(method);
 
    if (!Check())
    {
        MSGUser()->MSG_ERROR("Inputs check failed ! The fitting is not processed ! ");
        return;
    }
 
    if (method == ProfileFitter::FitMethod::ProfileLikelihood)
        FitPLH();
 
    if (method == ProfileFitter::FitMethod::AnalyticalChiSquare)
        FitACS();
 
    if (method == ProfileFitter::FitMethod::AnalyticalBoostedProfileLikelihood)
    {
        FitACS();
        FitPLH();
    }
 
    PrintResult();
}
 
void ProfileFitter::PrintResult()
{
    MSGUser()->MSG_INFO("// ------------------------ //");
    MSGUser()->MSG_INFO("// Parameter Fitting Result //");
    MSGUser()->MSG_INFO("// ------------------------ //");
    MSGUser()->MSG_INFO("Name             Type");
    for (auto &p : Par_Map)
    {
        if (p.second.isFixed)
            continue;
        if (p.second.Type == ParameterType::Interest)
            MSGUser()->MSG_INFO(p.second.Name.TString::Data(), " POI ", p.second.Value, "   ", p.second.Error);
    }
    for (auto &p : Par_Map)
    {
        if (p.second.isFixed)
            continue;
        if (p.second.Type == ParameterType::Normalization)
            MSGUser()->MSG_INFO(p.second.Name.TString::Data(), " Normalization ", p.second.Value, "   ", p.second.Error);
    }
    for (auto &p : Par_Map)
    {
        if (p.second.isFixed)
            continue;
        if (p.second.Type == ParameterType::Nuisance)
            MSGUser()->MSG_INFO(p.second.Name.TString::Data(), " Nuisance ", p.second.Value, "   ", p.second.Error);
    }
    for (auto &p : Par_Map)
    {
        if (p.second.isFixed)
            continue;
        if (p.second.Type == ParameterType::StatGamma)
            MSGUser()->MSG_INFO(p.second.Name.TString::Data(), " Gamma ", p.second.Value, "   ", p.second.Error);
    }
}
 
void ProfileFitter::PrintInfo()
{
    MSGUser()->MSG_INFO("// ---------- //");
    MSGUser()->MSG_INFO("// Basic Info //");
    MSGUser()->MSG_INFO("//----------- //");
    MSGUser()->MSG_INFO("Number of Observables:                   ", Obs_Map.size());
    MSGUser()->MSG_INFO("Number of Parameters:                    ", Par_Map.size());
    MSGUser()->MSG_INFO("Number of Parameters of Interest (POIS): ", GetNParameters(ParameterType::Interest));
    MSGUser()->MSG_INFO("Number of Nuisance Parameters:           ", GetNParameters(ParameterType::Nuisance));
    MSGUser()->MSG_INFO("Number of Normalization Parameter:       ", GetNParameters(ParameterType::Normalization));
    MSGUser()->MSG_INFO("// -------------- //");
    MSGUser()->MSG_INFO("// Parameter List //");
    MSGUser()->MSG_INFO("// -------------- //");
    MSGUser()->MSG_INFO("Name             Type");
    for (auto &p : Par_Map)
    {
        if (p.second.Type == ParameterType::Interest)
            MSGUser()->MSG_INFO(p.second.Name.TString::Data(), "      POI");
    }
    for (auto &p : Par_Map)
    {
        if (p.second.Type == ParameterType::Nuisance)
            MSGUser()->MSG_INFO(p.second.Name.TString::Data(), "      Nuisance");
    }
    for (auto &p : Par_Map)
    {
        if (p.second.Type == ParameterType::Normalization)
            MSGUser()->MSG_INFO(p.second.Name.TString::Data(), "      Normalization");
    }
    MSGUser()->MSG_INFO("// --------------- //");
    MSGUser()->MSG_INFO("// Observable List //");
    MSGUser()->MSG_INFO("// --------------- //");
    for (auto &obs : Obs_Map)
    {
        MSGUser()->MSG_INFO(obs.second.Name.TString::Data(), " NBins: ", obs.second.Data.size(), " NSamples: ", obs.second.Sample_Map.size());
        for (auto &samp : obs.second.Sample_Map)
            MSGUser()->MSG_INFO("     ", samp.second.Name, " NParameters: ", samp.second.Variation_Map.size(), " NVariations: ", samp.second.N_Variations);
    }
}
 
bool ProfileFitter::Check()
{
    // TODO
    return true;
}
 
void ProfileFitter::Prepare(ProfileFitter::FitMethod method)
{
    auto msgguard = MSGUser()->StartTitleWithGuard("ProfileFitter::Prepare");
 
    for (auto &obs : Obs_Map)
    {
        for (auto &samp : obs.second.Sample_Map)
        {
            samp.second.Start_Index = obs.second.Start_Index;
            samp.second.End_Index = obs.second.End_Index;
            samp.second.Rebin_Index = obs.second.Rebin_Index;
        }
    }
 
    for (auto &obs : Obs_Map)
        obs.second.SortVariations();
    if (method == FitMethod::ProfileLikelihood ||
        method == FitMethod::AnalyticalBoostedProfileLikelihood)
    {
        // Book Gamma Variations
        for (auto &obs : Obs_Map)
        {
            if (obs.second.isgamma_booked == true)
            {
                MSGUser()->MSG_ERROR("Statistical Nuisance Parameters for Observable ", obs.first.TString::Data(), " have already been booked !");
                return;
            }
 
            for (auto &g : Gamma_Map)
            {
                bool needtobook = false;
                vector<Parameter *> gamma_vec;
                if (g.second.size() > 0)
                {
                    for (int i = 0; i < g.second.size(); i++)
                    {
                        auto findsample = obs.second.Sample_Map.find(g.second[i]);
                        if (findsample != obs.second.Sample_Map.end())
                        {
                            needtobook = true;
                        }
                    }
                }
                if (needtobook)
                {
                    int ncount = 0;
                    Parameter *gamma_ptr;
                    for (int i = 0; i < obs.second.Data.size(); i++)
                    {
                        if ((i >= (obs.second.Start_Index > 0 ? obs.second.Start_Index : 0)) &&
                            (i <= (obs.second.End_Index > 0 ? obs.second.End_Index : obs.second.Data.size())))
                        {
                            if (ncount == 0)
                            {
                                TString gamma_name = TString::Format("Stat_%s_%s_Bin_%d", g.first.TString::Data(), obs.first.TString::Data(), i);
                                BookParameter(gamma_name, 0, -MAX_NUISANCE_SIGMA, MAX_NUISANCE_SIGMA, ParameterType::StatGamma);
                                MSGUser()->MSG_INFO("Book Stat Gamma: ", gamma_name.TString::Data());
                                gamma_ptr = GetParameter(gamma_name);
                                gamma_vec.emplace_back(gamma_ptr);
                                ncount = ncount + 1;
                            }
                            else
                            {
                                gamma_vec.emplace_back(gamma_ptr);
                                ncount = ncount + 1;
                            }
                            if (ncount == obs.second.Rebin_Index)
                                ncount = 0;
                        }
                        else
                            gamma_vec.emplace_back(nullptr);
                    }
 
                    for (auto &samp : obs.second.Sample_Map)
                    {
                        bool isprefix = false;
                        if (g.second.size() > 0)
                        {
                            for (int i = 0; i < g.second.size(); i++)
                            {
                                if (g.second[i].TString::EqualTo(samp.first.TString::Data()))
                                    isprefix = true;
                            }
                        }
                        if (isprefix)
                            samp.second.Gamma_Vec = gamma_vec;
                    }
                }
            }
            obs.second.isgamma_booked = true;
        }
 
        // Calculate the sigma of the gammas
        for (auto &par : Par_Map)
        {
            if (par.second.Type != ParameterType::StatGamma)
                continue;
 
            Parameter *par_ptr = GetParameter(par.first);
 
            double val = 0;
            double err = 0;
 
            for (auto &obs : Obs_Map)
            {
                for (auto &samp : obs.second.Sample_Map)
                {
                    for (int index = 0; index < samp.second.Gamma_Vec.size(); index++)
                    {
                        if (par_ptr == samp.second.Gamma_Vec[index])
                        {
                            val = val + samp.second.Nominal[index];
                            err = err + pow(samp.second.Nominal_Stat[index], 2);
                        }
                    }
                }
            }
 
            err = sqrt(err);
 
            par_ptr->Nuisance_Gaus_Mean = 1;
            par_ptr->Nuisance_Gaus_Sigma = err / val;
            MSGUser()->MSG_INFO(par_ptr->Name.TString::Data(),
                                " Mean : ", par_ptr->Nuisance_Gaus_Mean,
                                " Sigma : ", par_ptr->Nuisance_Gaus_Sigma);
 
            par_ptr->Init = 1;
            par_ptr->Value = 1;
            par_ptr->Error = err / val;
            par_ptr->Max = 1 + MAX_NUISANCE_SIGMA * err / val;
            par_ptr->Min = 1 - MAX_NUISANCE_SIGMA * err / val;
        }
    }
 
    // Link Normalization Factor
    for (auto &obs : Obs_Map)
    {
        for (auto &samp : obs.second.Sample_Map)
        {
            for (auto &par : Norm_Map)
            {
                if (par.second.size() > 0)
                {
                    for (int i = 0; i < par.second.size(); i++)
                    {
                        if (par.second[i].TString::EqualTo(samp.second.Name.TString::Data()))
                        {
                            samp.second.Norm_Par = GetParameter(par.first);
                        }
                    }
                }
            }
        }
    }
 
    // Symmetrization
    for (auto &obs : Obs_Map)
    {
        for (auto &samp : obs.second.Sample_Map)
        {
            for (auto &par : samp.second.NeedSymmetrization_Map)
            {
                if (par.second == true)
                {
                    MSGUser()->MSG_INFO("Symmetrize Variations for Obs: ", obs.first.TString::Data(), " Sample: ", samp.first.TString::Data(), " Parameter: ", par.first->Name.TString::Data());
                    auto findnom = samp.second.Variation_Nominal_Map.find(par.first);
                    auto findup = samp.second.Variation_Up_Map.find(par.first);
                    auto finddown = samp.second.Variation_Down_Map.find(par.first);
                    if (findnom != samp.second.Variation_Nominal_Map.end())
                    {
                        for (int i = 0; i < samp.second.Variation_Sigma_Map[par.first].size(); i++)
                        {
                            for (int j = 0; j < samp.second.Variation_Sigma_Map[par.first].size(); j++)
                            {
                                if (fabs(samp.second.Variation_Sigma_Map[par.first][i] + samp.second.Variation_Sigma_Map[par.first][j]) < 1e-5)
                                {
                                    if (samp.second.Variation_Sigma_Map[par.first][i] > 0)
                                        SymmetrizeVariation(samp.second.Variation_Map[par.first][i],
                                                            samp.second.Variation_Map[par.first][j],
                                                            samp.second.Variation_Nominal_Map[par.first]);
                                }
                            }
                        }
                        if ((findup != samp.second.Variation_Up_Map.end()) &&
                            (finddown != samp.second.Variation_Down_Map.end()))
                            SymmetrizeVariation(samp.second.Variation_Up_Map[par.first],
                                                samp.second.Variation_Down_Map[par.first],
                                                samp.second.Variation_Nominal_Map[par.first]);
                    }
                }
            }
        }
    }
 
    // Smooth
    for (auto &obs : Obs_Map)
    {
        for (auto &samp : obs.second.Sample_Map)
        {
            for (auto &par : samp.second.SmoothLevel_Map)
            {
                if (par.second > 0)
                {
                    MSGUser()->MSG_INFO("Smooth Variations for Obs: ", obs.first.TString::Data(), " Sample: ", samp.first.TString::Data(), " Parameter: ", par.first->Name.TString::Data());
                    auto findnom = samp.second.Variation_Nominal_Map.find(par.first);
                    auto findup = samp.second.Variation_Up_Map.find(par.first);
                    auto finddown = samp.second.Variation_Down_Map.find(par.first);
                    auto findvar = samp.second.Variation_Map.find(par.first);
                    if (findnom != samp.second.Variation_Nominal_Map.end())
                    {
                        if (findup != samp.second.Variation_Up_Map.end())
                            SmoothVariation(findup->second, findnom->second, par.second);
                        if (finddown != samp.second.Variation_Down_Map.end())
                            SmoothVariation(finddown->second, findnom->second, par.second);
 
                        if (findvar != samp.second.Variation_Map.end())
                        {
                            for (int i = 0; i < (findvar->second).size(); i++)
                                SmoothVariation((findvar->second)[i], findnom->second, par.second);
                        }
                    }
                }
            }
        }
    }
}
 
int ProfileFitter::GetNParameters(ParameterType type)
{
    int n = 0;
    if (type == ParameterType::Unknown)
        return Par_Map.size();
    for (auto &p : Par_Map)
    {
        if (p.second.Type == type)
            n = n + 1;
    }
    return n;
}
 
TString ProfileFitter::GetParameterName(int i)
{
    if (i >= Par_Map.size() || i < 0)
    {
        MSGUser()->MSG_ERROR("Negative Index or Index Too Large! Return Parameter Name : Unknown");
        return "Unknown";
    }
    else
    {
        int n = 0;
        for (auto &p : Par_Map)
        {
            if (n == i)
                return p.second.Name;
            n = n + 1;
        }
    }
    return "Unknown";
}
 
double ProfileFitter::Observable::Sample::GetProfileBinVariation(int i, double aim_sigma, Parameter *par)
{
    // if (std::isnan(aim_sigma))
    //     aim_sigma = 0.0;
 
    vector<double> sigma_vec = Variation_Sigma_Map[par];
    vector<vector<double>> variation_vec = Variation_Map[par];
    std::shared_ptr<ProfileCache> m_cache = Variation_Cache[par][i];
 
    if (sigma_vec.size() < 2 || variation_vec.size() != sigma_vec.size())
        cout << "Variation Number Less Than 2 !" << endl;
 
    if (par->Strategy == Parameter::ExtrapStrategy::POLY6EXP)
    {
 
        bool isok = true;
 
        auto findnom = Variation_Nominal_Map.find(par);
        auto finddown = Variation_Down_Map.find(par);
        auto findup = Variation_Up_Map.find(par);
        if (findnom != Variation_Nominal_Map.end() &&
            (finddown != Variation_Down_Map.end() || findup != Variation_Up_Map.end()))
            isok = true;
        else
        {
            cout << par->Name.TString::Data() << " Not Satisfy the condition of Extrap Strategy POLY6EXP !" << endl;
            isok = false;
        }
 
        // 1 + a1 x + a2 x^2 + a3 x^3 + ... a6 x^6
        // a1 + 2 a2 x + ... + 6 a6 x^5
        // 2 a2 + 3 * 2 a3 + ... + 6 * 5 a6 x^4
        if (isok)
        {
 
            if (m_cache == nullptr)
            {
                m_cache = std::make_shared<ProfileCache>();
                Variation_Cache[par][i] = m_cache;
 
                double base = 1;
                if (findup != Variation_Up_Map.end())
                    base = (findup->second)[i] / (findnom->second)[i];
                else
                    base = (-(finddown->second)[i] + 2 * (findnom->second)[i]) / (findnom->second)[i];
                m_cache->SaveCache("BaseUp", base);
 
                if (finddown != Variation_Down_Map.end())
                    base = (finddown->second)[i] / (findnom->second)[i];
                else
                    base = (-(findup->second)[i] + 2 * (findnom->second)[i]) / (findnom->second)[i];
                m_cache->SaveCache("BaseDown", base);
 
                double Y1 = 0;
                if (findup != Variation_Up_Map.end())
                    Y1 = (findup->second)[i];
                else
                    Y1 = (-(finddown->second)[i] + 2 * (findnom->second)[i]);
                double Y2 = 0;
                if (finddown != Variation_Down_Map.end())
                    Y2 = (finddown->second)[i];
                else
                    Y2 = (-(findup->second)[i] + 2 * (findnom->second)[i]);
                double Y0 = (findnom->second)[i];
 
                double LGY1 = TMath::Log(Y1 / Y0);
                double LGY2 = TMath::Log(Y2 / Y0);
                double LGY1S = pow(TMath::Log(Y1 / Y0), 2);
                double LGY2S = pow(TMath::Log(Y2 / Y0), 2);
 
                double a[6];
                a[0] = -(-15 * Y1 + 15 * Y2 + 7 * Y1 * LGY1 - Y1 * LGY1S - 7 * Y2 * LGY2 + Y2 * LGY2S) / (16 * Y0);
                a[1] = -(48 * Y0 - 24 * Y1 - 24 * Y2 + 9 * Y1 * LGY1 - Y1 * LGY1S + 9 * Y2 * LGY2 - Y2 * LGY2S) / (16 * Y0);
                a[2] = -(5 * Y1 - 5 * Y2 - 5 * Y1 * LGY1 + Y1 * LGY1S + 5 * Y2 * LGY2 - Y2 * LGY2S) / (8 * Y0);
                a[3] = -(-24 * Y0 + 12 * Y1 + 12 * Y2 - 7 * Y1 * LGY1 + Y1 * LGY1S - 7 * Y2 * LGY2 + Y2 * LGY2S) / (8 * Y0);
                a[4] = -(-3 * Y1 + 3 * Y2 + 3 * Y1 * LGY1 - Y1 * LGY1S - 3 * Y2 * LGY2 + Y2 * LGY2S) / (16 * Y0);
                a[5] = -(16 * Y0 - 8 * Y1 - 8 * Y2 + 5 * Y1 * LGY1 - Y1 * LGY1S + 5 * Y2 * LGY2 - Y2 * LGY2S) / (16 * Y0);
 
                m_cache->SaveCache("A0", a[0]);
                m_cache->SaveCache("A1", a[1]);
                m_cache->SaveCache("A2", a[2]);
                m_cache->SaveCache("A3", a[3]);
                m_cache->SaveCache("A4", a[4]);
                m_cache->SaveCache("A5", a[5]);
            }
 
            if (aim_sigma >= 1)
                return (findnom->second)[i] * pow(m_cache->GetCache("BaseUp"), aim_sigma);
            else if (aim_sigma <= -1)
                return (findnom->second)[i] * pow(m_cache->GetCache("BaseDown"), -aim_sigma);
            else
            {
 
                double a[6];
 
                a[0] = m_cache->GetCache("A0");
                a[1] = m_cache->GetCache("A1");
                a[2] = m_cache->GetCache("A2");
                a[3] = m_cache->GetCache("A3");
                a[4] = m_cache->GetCache("A4");
                a[5] = m_cache->GetCache("A5");
 
                double base = aim_sigma;
                double sum = 1.0;
                for (int k = 0; k < 6; k++)
                {
                    sum = sum + a[k] * base;
                    base = base * aim_sigma;
                }
 
                return (findnom->second)[i] * sum;
            }
        }
    }
 
    // Linear Extrapolation
 
    if (m_cache == nullptr)
    {
        m_cache = std::make_shared<ProfileCache>();
        Variation_Cache[par][i] = m_cache;
 
        for (int index = 0; index < sigma_vec.size() - 1; index++)
        {
            double X1 = sigma_vec[index];
            double X2 = sigma_vec[index + 1];
            double Y1 = variation_vec[index][i];
            double Y2 = variation_vec[index + 1][i];
 
            double A = Y2 - Y1;
            double B = X1 - X2;
            double C = X2 * Y1 - X1 * Y2;
 
            m_cache->SaveCache(TString::Format("A_%d", index), A);
            m_cache->SaveCache(TString::Format("B_%d", index), B);
            m_cache->SaveCache(TString::Format("C_%d", index), C);
        }
    }
 
    int index1 = -1;
    int index2 = -1;
    if (aim_sigma <= sigma_vec[0])
    {
        index1 = 0;
        index2 = 1;
    }
    else if (aim_sigma >= sigma_vec[sigma_vec.size() - 1])
    {
        index1 = sigma_vec.size() - 1;
        index2 = sigma_vec.size() - 2;
    }
    else
    {
        for (int j = 0; j < sigma_vec.size() - 1; j++)
        {
            if (aim_sigma >= sigma_vec[j] &&
                aim_sigma < sigma_vec[j + 1])
            {
                index1 = j;
                index2 = j + 1;
            }
        }
    }
 
    if (index1 < 0 || index2 < 0)
    {
        cout << "Negative Extrapolation Index for parameter " << par->Name << endl;
        std::cout << aim_sigma << std::endl;
        for (int i = 0; i < sigma_vec.size(); i++)
            std::cout << sigma_vec[i] << " ";
        std::cout << std::endl;
    }
 
    double A;
    double B;
    double C;
 
    if (index1 < index2)
    {
        A = m_cache->GetCache(TString::Format("A_%d", index1));
        B = m_cache->GetCache(TString::Format("B_%d", index1));
        C = m_cache->GetCache(TString::Format("C_%d", index1));
    }
    else
    {
        A = m_cache->GetCache(TString::Format("A_%d", index2));
        B = m_cache->GetCache(TString::Format("B_%d", index2));
        C = m_cache->GetCache(TString::Format("C_%d", index2));
    }
 
    return (-A * aim_sigma - C) / B;
}
 
void ProfileFitter::Observable::GetProfile(std::map<TString, Parameter> &m_par)
{
    for (auto &samp : Sample_Map)
        samp.second.GetProfile(m_par);
 
    for (int i = (Start_Index > 0 ? Start_Index : 0); i < (End_Index > 0 ? End_Index + 1 : Profile.size()); i++)
    {
        double sum = 0;
        for (auto &samp : Sample_Map)
            sum = sum + samp.second.Profile[i];
        Profile[i] = sum;
    }
}
 
void ProfileFitter::Observable::Sample::GetProfile(std::map<TString, Parameter> &m_par)
{
    for (int i = (Start_Index > 0 ? Start_Index : 0); i < (End_Index > 0 ? End_Index + 1 : Profile.size()); i++)
    {
        Profile[i] = Nominal[i];
 
        for (auto &itr : Variation_Sigma_Map)
        {
            if (Variation_Cache[itr.first][i] == nullptr)
            {
                double aim_profile = GetProfileBinVariation(i, (m_par[itr.first->Name]).Value, itr.first);
                double aim_nominal = GetProfileBinVariation(i, (m_par[itr.first->Name]).Init, itr.first);
                Variation_Cache[itr.first][i]->SaveCache("Nominal", aim_nominal);
                Profile[i] = Profile[i] + (aim_profile - aim_nominal);
 
                // std::cout << this->Name << " " << itr.first->Name << " " << i << " " << aim_profile << " " << aim_nominal << std::endl;
            }
            else
            {
                double aim_profile = GetProfileBinVariation(i, (m_par[itr.first->Name]).Value, itr.first);
                double aim_nominal = Variation_Cache[itr.first][i]->GetCache("Nominal");
                Profile[i] = Profile[i] + (aim_profile - aim_nominal);
 
                // std::cout << this->Name << " " << itr.first->Name << " " << i << " " << aim_profile << " " << aim_nominal << std::endl;
            }
        }
 
        if (Norm_Par != nullptr)
        {
            Profile[i] = Profile[i] * (m_par[Norm_Par->Name]).Value;
        }
 
        if (Gamma_Vec.size() > 0)
        {
            Profile[i] = Profile[i] * (m_par[Gamma_Vec[i]->Name]).Value;
        }
    }
}
 
double ProfileFitter::Observable::GetLikelihood()
{
    double prob = 0;
 
    int ncount = 0;
    double D = 0;
    double P = 0;
    double D_S = 0;
 
    for (int i = (Start_Index > 0 ? Start_Index : 0); i < (End_Index > 0 ? End_Index + 1 : Profile.size()); i++)
    {
        ncount = ncount + 1;
        D = D + Data[i];
        P = P + Profile[i];
        if (!USE_POISSON)
            D_S = D_S + Data_Stat[i] * Data_Stat[i];
 
        if (ncount == Rebin_Index)
        {
            double p = 1;
            if (USE_POISSON)
                p = D * std::log(P) - P - std::lgamma(D + 1.);
            else
                p = -0.5 * std::pow((D - P) / std::sqrt(D_S), 2);
 
            prob += p;
 
            ncount = 0;
            D = 0;
            P = 0;
            D_S = 0;
        }
    }
 
    if (ncount > 0)
    {
        double p = 1;
        if (USE_POISSON)
            p = D * std::log(P) - P - std::lgamma(D + 1.);
        else
            p = -0.5 * std::pow((D - P) / std::sqrt(D_S), 2);
 
        prob += p;
    }
 
    return prob;
}
 
void ProfileFitter::Observable::SortVariations()
{
    for (auto &samp : Sample_Map)
        samp.second.SortVariations();
}
 
void ProfileFitter::Observable::Sample::SortVariations()
{
    std::map<Parameter *, std::vector<double>>::iterator sig_itr;
    std::map<Parameter *, std::vector<std::vector<double>>>::iterator var_itr;
 
    var_itr = Variation_Map.begin();
    for (sig_itr = Variation_Sigma_Map.begin(); sig_itr != Variation_Sigma_Map.end(); sig_itr++)
    {
        if (sig_itr->second.size() >= 2)
        {
            for (int i = 0; i < sig_itr->second.size() - 1; i++)
            {
                for (int j = i + 1; j < sig_itr->second.size(); j++)
                {
                    if ((sig_itr->second)[i] > (sig_itr->second)[j])
                    {
                        double temp = (sig_itr->second)[i];
                        (sig_itr->second)[i] = (sig_itr->second)[j];
                        (sig_itr->second)[j] = temp;
 
                        std::vector<double> temp_vec = (var_itr->second)[i];
                        (var_itr->second)[i] = (var_itr->second)[j];
                        (var_itr->second)[j] = temp_vec;
                    }
                }
            }
        }
        var_itr++;
    }
}
 
void ProfileFitter::FitACS()
{
    // Calculate Gamma Matrix
    // Sensitivity Matrix : dt_i / da_l
    int Total_Template_BinNum = 0;
    for (auto &obs : Obs_Map)
    {
        Total_Template_BinNum = Total_Template_BinNum +
                                (obs.second.End_Index > 0 ? obs.second.End_Index + 1 : obs.second.GetNbins()) -
                                (obs.second.Start_Index > 0 ? obs.second.Start_Index : 0);
    }
    int NP_Num = 0;
    int POI_Num = 0;
    for (auto &par : Par_Map)
    {
        if (par.second.isFixed)
            continue;
 
        if (par.second.Type == ParameterType::Interest)
            POI_Num++;
        if (par.second.Type == ParameterType::Nuisance)
            NP_Num++;
        if (par.second.Type == ParameterType::Normalization)
            POI_Num++;
    }
 
    Gamma_Matrix = Eigen::MatrixXd(Total_Template_BinNum, NP_Num);
    V_Matrix = Eigen::MatrixXd(Total_Template_BinNum, Total_Template_BinNum);
    H_Matrix = Eigen::MatrixXd(Total_Template_BinNum, POI_Num);
    Y0_Matrix = Eigen::MatrixXd(Total_Template_BinNum, 1);
 
    int bin_index = 0;
    int par_index = 0;
    int poi_index = 0;
    for (auto &obs : Obs_Map)
    {
        for (int index = (obs.second.Start_Index > 0 ? obs.second.Start_Index : 0); index < (obs.second.End_Index > 0 ? obs.second.End_Index + 1 : obs.second.GetNbins()); index++)
        {
            par_index = 0;
            poi_index = 0;
            for (auto &par : Par_Map)
            {
                if (par.second.isFixed)
                    continue;
 
                if ((par.second.Type == ParameterType::Nuisance) ||
                    (par.second.Type == ParameterType::Interest))
                {
                    Parameter *par_ptr = &(par.second);
                    double sum_k = 0;
                    for (auto &samp : obs.second.Sample_Map)
                    {
                        if (samp.second.Variation_Map.find(par_ptr) != samp.second.Variation_Map.end())
                        {
                            double k = 0;
                            // calculate k
 
                            // Average values
                            double ax = 0;
                            double ay = 0;
                            for (int v = 0; v < samp.second.Variation_Sigma_Map[par_ptr].size(); v++)
                            {
                                ax = ax + samp.second.Variation_Sigma_Map[par_ptr][v];
                                ay = ay + samp.second.Variation_Map[par_ptr][v][index];
                            }
                            ax = ax / samp.second.Variation_Sigma_Map[par_ptr].size();
                            ay = ay / samp.second.Variation_Sigma_Map[par_ptr].size();
 
                            double top = 0;
                            double bottom = 0;
                            for (int v = 0; v < samp.second.Variation_Sigma_Map[par_ptr].size(); v++)
                            {
                                top = top + (samp.second.Variation_Sigma_Map[par_ptr][v] - ax) * (samp.second.Variation_Map[par_ptr][v][index] - ay);
                                bottom = bottom + pow(samp.second.Variation_Sigma_Map[par_ptr][v] - ax, 2);
                            }
 
                            k = top / bottom;
 
                            sum_k = sum_k + k;
                        }
                    }
 
                    if (!isfinite(sum_k))
                        MSGUser()->MSG_INFO("Sensitivity : ", sum_k);
                    if (par.second.Type == ParameterType::Nuisance)
                    {
                        Gamma_Matrix(bin_index, par_index) = sum_k;
                        par_index++;
                    }
                    if (par.second.Type == ParameterType::Interest)
                    {
                        H_Matrix(bin_index, poi_index) = sum_k;
                        poi_index++;
                    }
                }
 
                if (par.second.Type == ParameterType::Normalization)
                {
                    Parameter *par_ptr = &(par.second);
                    double sum_k = 0;
                    for (auto &samp : obs.second.Sample_Map)
                    {
                        if (samp.second.Norm_Par == par_ptr)
                        {
                            sum_k = sum_k + samp.second.Nominal[index];
                        }
                    }
 
                    if (!isfinite(sum_k))
                        MSGUser()->MSG_INFO("Sensitivity : ", sum_k);
 
                    H_Matrix(bin_index, poi_index) = sum_k;
                    poi_index++;
                }
            }
 
            for (int i = 0; i < Total_Template_BinNum; i++)
            {
                V_Matrix(bin_index, i) = 0;
                V_Matrix(i, bin_index) = 0;
            }
            double stat_err = pow(obs.second.Data_Stat[index], 2);
            double y0 = obs.second.Data[index];
            for (auto &samp : obs.second.Sample_Map)
            {
                stat_err = stat_err + pow(samp.second.Nominal_Stat[index], 2);
                y0 = y0 - samp.second.Nominal[index];
            }
            V_Matrix(bin_index, bin_index) = stat_err;
            Y0_Matrix(bin_index, 0) = y0;
            if (!isfinite(stat_err))
                MSGUser()->MSG_INFO("Stat_Err^2 : ", stat_err);
 
            bin_index++;
        }
    }
 
    Q_Matrix = (Eigen::MatrixXd::Identity(NP_Num, NP_Num) + Gamma_Matrix.transpose() * V_Matrix.inverse() * Gamma_Matrix).inverse() * (Gamma_Matrix.transpose() * V_Matrix.inverse());
    S_Matrix = V_Matrix.inverse() * (Eigen::MatrixXd::Identity(Total_Template_BinNum, Total_Template_BinNum) - Gamma_Matrix * Q_Matrix);
    C_Matrix = V_Matrix + Gamma_Matrix * Gamma_Matrix.transpose();
    Lambda_Matrix = (-1) * (H_Matrix.transpose() * S_Matrix * H_Matrix).inverse() * H_Matrix.transpose() * S_Matrix;
    POI_Shift = Lambda_Matrix * Y0_Matrix;
    Nuisance_Shift = Q_Matrix * (Eigen::MatrixXd::Identity(Total_Template_BinNum, Total_Template_BinNum) + H_Matrix * Lambda_Matrix) * Y0_Matrix;
 
    POI_Cov_Matrix = (H_Matrix.transpose() * S_Matrix * H_Matrix).inverse();
 
    Rho_Matrix = (H_Matrix.transpose() * V_Matrix.inverse() * H_Matrix).inverse() * (H_Matrix.transpose() * V_Matrix.inverse());
    Epsilon_Matrix = H_Matrix * Rho_Matrix - Eigen::MatrixXd::Identity(Total_Template_BinNum, Total_Template_BinNum);
    Nuisance_Cov_Matrix = (Eigen::MatrixXd::Identity(NP_Num, NP_Num) + (Epsilon_Matrix * Gamma_Matrix).transpose() * V_Matrix.inverse() * (Epsilon_Matrix * Gamma_Matrix)).inverse();
 
    poi_index = 0;
    par_index = 0;
    for (auto &par : Par_Map)
    {
        if (par.second.isFixed)
            continue;
        if (par.second.Type == ParameterType::Nuisance)
        {
            par.second.Value = Nuisance_Shift(par_index, 0);
            par.second.Error = sqrt(Nuisance_Cov_Matrix(par_index, par_index));
            par_index++;
        }
        if (par.second.Type == ParameterType::Normalization)
        {
            par.second.Value = par.second.Init - POI_Shift(poi_index, 0);
            par.second.Error = sqrt(POI_Cov_Matrix(poi_index, poi_index));
            poi_index++;
        }
        if (par.second.Type == ParameterType::Interest)
        {
            par.second.Value = par.second.Init - POI_Shift(poi_index, 0);
            par.second.Error = sqrt(POI_Cov_Matrix(poi_index, poi_index));
            poi_index++;
        }
    }
}
 
void ProfileFitter::FitPLH()
{
    ProfileFitter::PLH_FCN fcn(pool, MSGUser(), &Par_Map, &Obs_Map);
 
    // Define Parameters
    ROOT::Minuit2::MnUserParameters upar;
    for (auto &p : Par_Map)
    {
        if (p.second.Type == ParameterType::Interest)
        {
            upar.Add(p.second.Name.TString::Data(),
                     p.second.Value,
                     15,
                     p.second.Min,
                     p.second.Max);
        }
        else
        {
            upar.Add(p.second.Name.TString::Data(),
                     p.second.Value,
                     1,
                     p.second.Min,
                     p.second.Max);
        }
        if (p.second.isFixed)
            upar.Fix(p.second.Name.TString::Data());
    }
 
    // Pre Run
    int index = 0;
    vector<int> NuisanceIndex_vec;
    vector<double> DeltaLH_vec;
    for (auto &p : Par_Map)
    {
        if ((p.second.Type == ParameterType::Nuisance) &&
            (!p.second.isFixed))
        {
            ROOT::Minuit2::MnScan scan(fcn, upar, 2);
            auto res = scan.Scan(index, 3, -1, 1);
            double delta_lh = res[1].second - res[2].second;
            if (fabs(delta_lh - 0.5) < 1e-5)
            {
                upar.Fix(index);
                MSGUser()->MSG_INFO(p.second.Name.TString::Data(), " Fixed to 0 ! ( ", fabs(delta_lh - 0.5), " )");
            }
            else
            {
                NuisanceIndex_vec.emplace_back(index);
                DeltaLH_vec.emplace_back(fabs(delta_lh - 0.5));
            }
        }
        index = index + 1;
    }
 
    if (DeltaLH_vec.size() > 1)
    {
        for (int i = 0; i < DeltaLH_vec.size() - 1; i++)
        {
            for (int j = 1; j < DeltaLH_vec.size(); j++)
            {
                if (DeltaLH_vec[i] < DeltaLH_vec[j])
                {
                    int temp_index = NuisanceIndex_vec[i];
                    NuisanceIndex_vec[i] = NuisanceIndex_vec[j];
                    NuisanceIndex_vec[j] = temp_index;
                    double temp_delta = DeltaLH_vec[i];
                    DeltaLH_vec[i] = DeltaLH_vec[j];
                    DeltaLH_vec[j] = temp_delta;
                }
            }
        }
    }
 
    MSGUser()->MSG_INFO("Number of Nuisance Parameters to Fit: ", NuisanceIndex_vec.size());
    // for (int i = 20; i < NuisanceIndex_vec.size(); i++)
    //     upar.Fix(NuisanceIndex_vec[i]);
 
    fcn.Reset();
    MSGUser()->MSG_INFO("Start Simplex");
    ROOT::Minuit2::MnSimplex simplex(fcn, upar, 2);
    auto minsimplex = simplex();
    MSGUser()->MSG_INFO("Start Migrad");
    ROOT::Minuit2::MnMigrad migrad(fcn, minsimplex.UserState(), ROOT::Minuit2::MnStrategy(2));
    auto minmigrad = migrad();
 
    MSGUser()->MSG_INFO("Final -Log(LH) : ", minmigrad.Fval());
 
    index = 0;
    for (auto &p : Par_Map)
    {
        p.second.Value = minmigrad.UserState().Value(index);
        p.second.Error = minmigrad.UserState().Error(index);
        p.second.isFitted = true;
        index = index + 1;
    }
}
 
double ProfileFitter::PLH_FCN::operator()(const std::vector<double> &vpara) const
{
    std::map<TString, Parameter> Temp_Par_Map = (*Par_Map);
    int index = 0;
    double prob = 0;
 
    if (pool == nullptr)
    {
        index = 0;
        std::map<TString, Parameter>::iterator par_itr;
        for (par_itr = Temp_Par_Map.begin(); par_itr != Temp_Par_Map.end(); par_itr++)
        {
            par_itr->second.Value = vpara[index];
            prob = prob + par_itr->second.GetLikelihood();
            index = index + 1;
        }
 
        std::map<TString, Observable>::iterator obs_itr;
        for (obs_itr = Obs_Map->begin(); obs_itr != Obs_Map->end(); obs_itr++)
        {
            obs_itr->second.GetProfile(Temp_Par_Map);
            prob = prob + obs_itr->second.GetLikelihood();
        }
    } // Single Thread Mode
    else
    {
        index = 0;
        std::map<TString, Parameter>::iterator par_itr;
        for (par_itr = Temp_Par_Map.begin(); par_itr != Temp_Par_Map.end(); par_itr++)
        {
            par_itr->second.Value = vpara[index];
            index = index + 1;
        }
 
        std::vector<std::future<StatusCode>> vfuture_obs;
        std::map<TString, Observable>::iterator obs_itr;
        for (obs_itr = Obs_Map->begin(); obs_itr != Obs_Map->end(); obs_itr++)
        {
            for (auto &samp : obs_itr->second.Sample_Map)
            {
                vfuture_obs.emplace_back(pool->enqueue(Observable::Sample::GetProfile_Thread, &(samp.second), Temp_Par_Map));
            }
        }
 
        std::vector<std::future<double>> vfuture_par;
        for (auto &par : Temp_Par_Map)
            vfuture_par.emplace_back(pool->enqueue(Parameter::GetLikelihood_Thread, &(par.second)));
 
        for (int i = 0; i < vfuture_obs.size(); i++)
            vfuture_obs[i].get();
        vfuture_obs.clear();
 
        for (obs_itr = Obs_Map->begin(); obs_itr != Obs_Map->end(); obs_itr++)
            vfuture_par.emplace_back(pool->enqueue(Observable::GetLikelihood_Thread, &(obs_itr->second)));
 
        for (int i = 0; i < vfuture_par.size(); i++)
            prob = prob + vfuture_par[i].get();
        vfuture_par.clear();
 
    } // Multi Thread Mode
 
    if (NCall == 0)
    {
        MSGUser()->MSG_DEBUG("First -Log(LH): ", -prob);
        Max_Log_Likelihood = -prob;
    }
    else if ((-prob) < Max_Log_Likelihood)
    {
        PrintMutex.lock();
        index = 0;
        std::map<TString, Parameter>::iterator par_itr;
        for (par_itr = Temp_Par_Map.begin(); par_itr != Temp_Par_Map.end(); par_itr++)
        {
            if (par_itr->second.Type == ParameterType::Interest)
                MSGUser()->MSG_DEBUG("POI: ", par_itr->second.Name.TString::Data(), " = ", par_itr->second.Value, " ( Index No. ", index, " )");
            else if (par_itr->second.Value != par_itr->second.Init)
                MSGUser()->MSG_DEBUG("Nuisance: ", par_itr->second.Name.TString::Data(), " = ", par_itr->second.Value, " ( Index No. ", index, " )");
            index = index + 1;
        }
        MSGUser()->MSG_DEBUG("New Minimum -Log(LH): ", -prob, " ( NCall = ", NCall + 1, " )");
        PrintMutex.unlock();
        Max_Log_Likelihood = -prob;
    }
    // else if ((NCall + 1) % 1000 == 0)
    else
        MSGUser()->MSG_DEBUG("New -Log(LH): ", -prob, " ( NCall = ", NCall + 1, " )");
    NCall = NCall + 1;
 
    return -prob;
}
 
double ProfileFitter::Parameter::GetLikelihood()
{
    if (Type == ParameterType::Nuisance ||
        Type == ParameterType::StatGamma)
        return -0.5 * std::pow((Value - Nuisance_Gaus_Mean) / Nuisance_Gaus_Sigma, 2);
    // return -0.5 * std::pow((Value - Nuisance_Gaus_Mean) / Nuisance_Gaus_Sigma, 2) - std::log(2.50662827463100024 * Nuisance_Gaus_Sigma);
 
    return 0;
}
 
void ProfileFitter::SetGammaPrefix(TString samplename, TString prefix)
{
    auto findprefix = Gamma_Map.find(prefix);
    if (findprefix != Gamma_Map.end())
    {
        bool findsample = false;
        if (findprefix->second.size() > 0)
        {
            for (int i = 0; i < findprefix->second.size(); i++)
            {
                if ((findprefix->second)[i].TString::EqualTo(samplename))
                    findsample = true;
            }
        }
        if (!findsample)
        {
            findprefix->second.emplace_back(samplename);
        }
    }
    else
    {
        vector<TString> Sample_Vec;
        Sample_Vec.emplace_back(samplename);
        for (auto &g : Gamma_Map)
        {
            if (g.second.size() > 0)
            {
                for (int i = 0; i < g.second.size(); i++)
                {
                    if ((g.second)[i].TString::EqualTo(samplename.TString::Data()))
                    {
                        g.second.erase(g.second.begin() + i);
                        i--;
                    }
                }
            }
        }
        Gamma_Map.emplace(std::pair<TString, vector<TString>>(prefix, Sample_Vec));
    }
}
 
void ProfileFitter::LinkNormalizationSample(TString name, TString samplename)
{
    bool exist = false;
    for (auto &par : Norm_Map)
    {
        if (par.second.size() > 0)
        {
            for (int i = 0; i < par.second.size(); i++)
            {
                if (par.second[i].TString::EqualTo(samplename))
                    exist = true;
            }
        }
    }
    if (exist)
    {
        MSGUser()->MSG_ERROR("Normalization Factor for ", samplename.TString::Data(), " already exist!");
        return;
    }
    if (GetParameter(name) != nullptr)
    {
        auto findpar = Norm_Map.find(name);
        if (findpar != Norm_Map.end())
        {
            bool findsamp = false;
            if (findpar->second.size() > 0)
            {
                for (int i = 0; i < findpar->second.size(); i++)
                {
                    if ((findpar->second)[i].TString::EqualTo(samplename))
                        findsamp = true;
                }
            }
            if (!findsamp)
                findpar->second.emplace_back(samplename);
        }
        else
        {
            vector<TString> sample_vec;
            sample_vec.emplace_back(samplename);
            Norm_Map.emplace(std::pair<TString, vector<TString>>(name, sample_vec));
        }
    }
}
void ProfileFitter::SymmetrizeVariation(vector<double> &up, vector<double> &down, vector<double> &norm)
{
    for (int i = 0; i < norm.size(); i++)
    {
        double delta_up = up[i] - norm[i];
        double delta_down = down[i] - norm[i];
 
        double new_delta_up = (delta_up - delta_down) / 2;
        double new_delta_down = -new_delta_up;
 
        up[i] = norm[i] + new_delta_up;
        down[i] = norm[i] + new_delta_down;
 
        if (up[i] < 0)
            up[i] = 1e-5;
        if (down[i] < 0)
            down[i] = 1e-5;
    }
}
 
void ProfileFitter::SmoothVariation(vector<double> &var, vector<double> &norm, int level)
{
    if (level <= 0)
        return;
    if (norm.size() > 2)
    {
        double sum_init = 0;
        sum_init = accumulate(var.begin(), var.end(), 0);
        vector<double> new_ratio_vec;
        for (int i = 1; i < norm.size() - 1; i++)
        {
            double ratio_p = var[i - 1] / norm[i - 1];
            double ratio_n = var[i + 1] / norm[i + 1];
 
            double ratio = var[i] / norm[i];
 
            if (!isfinite(ratio_p))
                ratio_p = 0;
            if (!isfinite(ratio_n))
                ratio_n = 0;
            if (!isfinite(ratio))
                ratio = 0;
 
            double new_ratio = (level * ratio + ratio_p + ratio_n) / (level + 2);
 
            new_ratio_vec.emplace_back(new_ratio);
        }
 
        for (int i = 1; i < norm.size() - 1; i++)
            var[i] = norm[i] * new_ratio_vec[i - 1];
 
        double sum_new = 0;
        sum_new = accumulate(var.begin(), var.end(), 0);
 
        double scale = sum_init / sum_new;
        for (int i = 0; i < norm.size(); i++)
            var[i] = var[i] * scale;
    }
}
 
void ProfileFitter::RebinObservable(TString name, int index)
{
    if (index < 2)
    {
        MSGUser()->MSG_ERROR("Rebin Index ", index, " is smaller than 2 !");
        return;
    }
 
    auto findobs = Obs_Map.find(name);
    if (findobs == Obs_Map.end())
    {
        MSGUser()->MSG_ERROR("Observable ", name.TString::Data(), " has not been defined !");
        return;
    }
 
    findobs->second.Rebin_Index = index;
}
 
void ProfileFitter::SetObservableRange(TString name, int start, int end)
{
    if (end < start)
    {
        MSGUser()->MSG_ERROR("Start Index ", start, "is larger than End Index ", end, " !");
        return;
    }
 
    auto findobs = Obs_Map.find(name);
    if (findobs == Obs_Map.end())
    {
        MSGUser()->MSG_ERROR("Observable ", name.TString::Data(), " has not been defined !");
        return;
    }
 
    findobs->second.Start_Index = start;
    findobs->second.End_Index = end;
}
 
void ProfileFitter::SetSymmetrization(TString obsname, TString samplename, TString parname, bool need)
{
    auto msgguard = MSGUser()->StartTitleWithGuard("ProfileFitter::SetSymmetrization");
    auto findobs = Obs_Map.find(obsname);
    if (findobs != Obs_Map.end())
    {
        auto findsample = findobs->second.Sample_Map.find(samplename);
        if (findsample != findobs->second.Sample_Map.end())
        {
            auto findpar = findsample->second.NeedSymmetrization_Map.find(GetParameter(parname));
            if (findpar != findsample->second.NeedSymmetrization_Map.end())
            {
                findpar->second = need;
            }
            else
                MSGUser()->MSG_ERROR("Sample ", samplename.TString::Data(),
                                     " of Observable ", obsname.TString::Data(),
                                     " Not Affected by Parameter ", parname.TString::Data(), " !");
        }
        else
            MSGUser()->MSG_ERROR("Sample ", samplename.TString::Data(),
                                 " of Observable ", obsname.TString::Data(), " Not Defined !");
    }
    else
        MSGUser()->MSG_ERROR("Observable ", obsname.TString::Data(), " Not Defined !");
}
 
void ProfileFitter::SetSmooth(TString obsname, TString samplename, TString parname, int slevel)
{
    auto msgguard = MSGUser()->StartTitleWithGuard("ProfileFitter::SetSymmetrization");
    auto findobs = Obs_Map.find(obsname);
    if (findobs != Obs_Map.end())
    {
        auto findsample = findobs->second.Sample_Map.find(samplename);
        if (findsample != findobs->second.Sample_Map.end())
        {
            auto findpar = findsample->second.SmoothLevel_Map.find(GetParameter(parname));
            if (findpar != findsample->second.SmoothLevel_Map.end())
            {
                findpar->second = slevel;
            }
            else
                MSGUser()->MSG_ERROR("Sample ", samplename.TString::Data(),
                                     " of Observable ", obsname.TString::Data(),
                                     " Not Affected by Parameter ", parname.TString::Data(), " !");
        }
        else
            MSGUser()->MSG_ERROR("Sample ", samplename.TString::Data(),
                                 " of Observable ", obsname.TString::Data(), " Not Defined !");
    }
    else
        MSGUser()->MSG_ERROR("Observable ", obsname.TString::Data(), " Not Defined !");
}
 
void ProfileFitter::SetParameterGroup(TString gname, TString parname)
{
    auto findgroup = ParGroup_Map.find(gname);
    if (findgroup != ParGroup_Map.end())
    {
        if (GetParameter(parname) != nullptr)
        {
            findgroup->second.emplace_back(parname);
        }
    }
    else
    {
        if (GetParameter(parname) != nullptr)
        {
            vector<TString> empty_vec;
            empty_vec.emplace_back(parname);
            ParGroup_Map.emplace(std::pair<TString, vector<TString>>(gname, empty_vec));
        }
    }
}
 
std::map<TString, double> ProfileFitter::EstimateUnc(TString g_name, FitMethod method, int toy_num)
{
    auto findgroup = ParGroup_Map.find(g_name);
    if (findgroup != ParGroup_Map.end())
    {
        return EstimateUnc(findgroup->second, method, toy_num);
    }
    else
    {
        auto guard = MSGUser()->StartTitleWithGuard("ProfileFitter::EstimateUnc");
        MSGUser()->MSG_ERROR("Parameter group ", g_name.TString::Data(), " is not defined, return an empty map!");
        return std::map<TString, double>();
    }
}
 
std::map<TString, double> ProfileFitter::EstimateUnc(const std::vector<TString> &par_vec, FitMethod method, int toy_num)
{
    std::map<TString, vector<double>> Original_Data_Map;
    for (auto &obs : Obs_Map)
    {
        Original_Data_Map.emplace(std::pair<TString, vector<double>>(obs.first, obs.second.Data));
    }
 
    TRandom3 *my_rnd = new TRandom3(12345);
    std::map<TString, double> Sum_Par_Map;
    std::map<TString, double> Sum_Par2_Map;
 
    for (auto &par : Par_Map)
    {
        Sum_Par_Map.emplace(std::pair<TString, double>(par.first, 0));
        Sum_Par2_Map.emplace(std::pair<TString, double>(par.first, 0));
    }
 
    for (int i = 0; i < toy_num; i++)
    {
        if (i % 100 == 0 && i != 0)
            MSGUser()->MSG_DEBUG("Produced ", i, " Toys ! ( Total Number = ", toy_num, " )");
        std::map<TString, Parameter> Temp_Par_Map = Par_Map;
        for (auto &par : Temp_Par_Map)
        {
            par.second.Value = par.second.Init;
        }
 
        for (int j = 0; j < par_vec.size(); j++)
        {
            Parameter *par = GetParameter(par_vec[j]);
            Temp_Par_Map[par_vec[j]].Value = my_rnd->Gaus(par->Nuisance_Gaus_Mean, par->Nuisance_Gaus_Sigma);
        } // Random Nuisance
 
        for (auto &obs : Obs_Map)
        {
            obs.second.GetProfile(Temp_Par_Map);
            for (int j = 0; j < obs.second.Data.size(); j++)
            {
                obs.second.Data[j] = Original_Data_Map[obs.first][j] - (obs.second.Profile[j]);
                for (auto &samp : obs.second.Sample_Map)
                {
                    obs.second.Data[j] = obs.second.Data[j] + samp.second.Nominal[j];
                }
            }
        } // Shifted Data
 
        if (method == ProfileFitter::FitMethod::ProfileLikelihood)
            FitPLH();
 
        if (method == ProfileFitter::FitMethod::AnalyticalChiSquare)
            FitACS();
 
        if (method == ProfileFitter::FitMethod::AnalyticalBoostedProfileLikelihood)
        {
            FitACS();
            FitPLH();
        }
 
        for (auto &par : Par_Map)
        {
            Sum_Par_Map[par.first] = Sum_Par_Map[par.first] + par.second.Value;
            Sum_Par2_Map[par.first] = Sum_Par2_Map[par.first] + pow(par.second.Value, 2);
        }
    }
    std::map<TString, double> Unc_Map;
 
    for (auto &par : Par_Map)
    {
        // add fabs to prevent numerical errors
        double unc = sqrt(std::fabs(Sum_Par2_Map[par.first] / toy_num - pow(Sum_Par_Map[par.first] / toy_num, 2)));
        Unc_Map.emplace(std::pair<TString, double>(par.first, unc));
    }
 
    for (auto &obs : Obs_Map)
    {
        obs.second.Data = Original_Data_Map[obs.first];
    }
 
    delete my_rnd;
 
    return Unc_Map;
}
 
double ProfileFitter::GetChi2()
{
    std::map<TString, Parameter> Temp_Par_Map = Par_Map;
 
    double Chi2 = 0;
    for (auto &obs : Obs_Map)
    {
        obs.second.GetProfile(Temp_Par_Map);
        for (int i = 0; i < obs.second.Profile.size(); i++)
        {
            double err2 = pow(obs.second.Data_Stat[i], 2);
            double diff2 = pow(obs.second.Data[i] - obs.second.Profile[i], 2);
            for (auto &samp : obs.second.Sample_Map)
            {
                if (samp.second.Gamma_Vec.size() == 0)
                    err2 = err2 + pow(samp.second.Nominal_Stat[i], 2);
            }
            Chi2 = Chi2 + diff2 / err2;
        }
    }
 
    return Chi2;
}
 
double ProfileFitter::GetChi2(const std::map<TString, double> &par_val)
{
    std::map<TString, Parameter> Temp_Par_Map = Par_Map;
 
    for (auto &val : par_val)
    {
        Parameter *par = GetParameter(val.first);
        if (par != nullptr)
            Temp_Par_Map[val.first].Value = val.second;
    }
 
    double Chi2 = 0;
    for (auto &obs : Obs_Map)
    {
        obs.second.GetProfile(Temp_Par_Map);
        for (int i = 0; i < obs.second.Profile.size(); i++)
        {
            double err2 = pow(obs.second.Data_Stat[i], 2);
            double diff2 = pow(obs.second.Data[i] - obs.second.Profile[i], 2);
            for (auto &samp : obs.second.Sample_Map)
            {
                if (samp.second.Gamma_Vec.size() == 0)
                    err2 = err2 + pow(samp.second.Nominal_Stat[i], 2);
            }
            Chi2 = Chi2 + diff2 / err2;
        }
    }
 
    return Chi2;
}