GenericSFApplier.h

Go to the documentation of this file.

#pragma once
 
#include "Core/Objects/interface/Format.h"
#include "Core/Objects/interface/Weight.h"
 
#include <correction.h>
#include <darwin.h>
 
#include <TFile.h>
#include <TH1.h>
 
#include <filesystem>
#include <map>
#include <utility>
#include <vector>
 
namespace DAS {
 
[[deprecated]]
inline std::unique_ptr<TH1> getHistSafe (std::unique_ptr<TFile>& f,
                                         const std::string& name)
{
    auto h = std::unique_ptr<TH1>(f->Get<TH1>(name.c_str()));
    if (!h) {
        std::cerr << name << std::endl;
        BOOST_THROW_EXCEPTION( Darwin::Exceptions::BadInput("Can't find required histogram", f) );
    }
    h->SetDirectory(nullptr);
    return h;
}
 
template<class Object, class... Context>
class GenericHistSFApplier
{
    std::unique_ptr<TFile> m_file; 
    bool m_correction,    
         m_uncertainties; 
 
    std::unique_ptr<TH1> m_nominal;
    std::map<std::string, std::unique_ptr<TH1>> m_bin_wise_uncs;
    std::map<std::string, std::unique_ptr<TH1>> m_global_uncs;
 
public:
    GenericHistSFApplier (const std::filesystem::path& filePath,
                          bool correction, bool uncertainties);
    virtual ~GenericHistSFApplier () noexcept = default;
 
    void operator() (Object& object, const Context&... ctx) const;
    void operator() (std::vector<Object>& objects, const Context&... ctx) const;
    std::vector<std::string> weightNames () const;
 
protected:
    enum Interpretation {
        UseBinContent,
        UseBinError,
    };
 
    void loadNominal (const std::string& histPath);
    void loadBinWiseUnc (const std::string& name,
                         const std::string& histPath,
                         Interpretation intp = UseBinContent);
    void loadGlobalUnc (const std::string& name,
                        const std::string& histPath,
                        Interpretation intp = UseBinContent);
 
    virtual bool passes (const Object& obj, const Context&... ctx) const = 0;
 
    virtual int binIndex (const Object& obj,
                          const Context&... ctx,
                          const std::unique_ptr<TH1>& hist) const = 0;
};
 
template<class Object, class... Context>
GenericHistSFApplier<Object, Context...>::GenericHistSFApplier (
        const std::filesystem::path& filePath,
        bool correction,
        bool uncertainties)
    : m_correction(correction)
    , m_uncertainties(uncertainties)
{
    namespace fs = std::filesystem;
 
    if (!correction)
        return;
 
    // Check that the file exists
    if (!fs::exists(filePath))
        BOOST_THROW_EXCEPTION(
            fs::filesystem_error("Not found",
                                 filePath,
                                 make_error_code(std::errc::no_such_file_or_directory)));
 
    m_file = std::make_unique<TFile>(filePath.c_str(), "READ");
}
 
 
template<class Object, class... Context>
void GenericHistSFApplier<Object, Context...>::operator() (
        Object& obj, const Context&... ctx) const
{
    auto& weights = obj.weights;
 
    // For safety
    if (weights.empty()) weights.emplace_back(DAS::Weight{1, 0});
 
    const float old = weights.front().v;
 
    // Apply the cut
    if (!passes(obj, ctx...)) {
        weights *= 0;
        return; // Skip applying weights
    }
 
    // Nominal correction
    if (!m_correction) return;
 
    const int bin = binIndex(obj, ctx..., m_nominal);
    const float sf = m_nominal->GetBinContent(bin);
    weights *= sf;
 
    // Systematics
    if (!m_uncertainties) return;
 
    // The order of the next two blocks must match weightNames()
    for (const auto& [_, hist] : m_bin_wise_uncs) {
        const float variation = hist->GetBinContent(bin);
        weights.push_back({old * (sf + variation), bin});
        weights.push_back({old * (sf - variation), bin});
    }
 
    for (const auto& [_, hist] : m_global_uncs) {
        const float variation = hist->GetBinContent(bin);
        weights.push_back({old * (sf + variation), 0});
        weights.push_back({old * (sf - variation), 0});
    }
}
 
template<class Object, class... Context>
void GenericHistSFApplier<Object, Context...>::operator() (
        std::vector<Object>& objects, const Context&... ctx) const
{
    for (auto& obj: objects) {
        (*this)(obj, ctx...);
    }
}
 
template<class Object, class... Context>
std::vector<std::string>
GenericHistSFApplier<Object, Context...>::weightNames () const
{
    std::vector<std::string> names;
    // The order of the next two blocks must match the call operator
    for (const auto& [name, _] : m_bin_wise_uncs) {
        names.push_back(name + SysUp);
        names.push_back(name + SysDown);
    }
    for (const auto& [name, _] : m_global_uncs) {
        names.push_back(name + SysUp);
        names.push_back(name + SysDown);
    }
    return names;
}
 
template<class Object, class... Context>
void GenericHistSFApplier<Object, Context...>::loadNominal (
        const std::string& histPath)
{
    if (!m_correction) return;
    m_nominal = getHistSafe(m_file, histPath);
}
 
template<class Object, class... Context>
void GenericHistSFApplier<Object, Context...>::loadBinWiseUnc (
        const std::string& name, const std::string& histPath, Interpretation intp)
{
    if (!m_correction || !m_uncertainties) return;
    auto hist = getHistSafe(m_file, histPath);
    if (intp == UseBinError) // We always store variations in the bin contents
        for (int i = 0; i <= hist->GetNcells(); ++i)
            hist->SetBinContent(i, hist->GetBinError(i));
    m_bin_wise_uncs[name] = std::move(hist);
}
 
template<class Object, class... Context>
void GenericHistSFApplier<Object, Context...>::loadGlobalUnc (
        const std::string& name, const std::string& histPath, Interpretation intp)
{
    if (!m_correction || !m_uncertainties) return;
    auto hist = getHistSafe(m_file, histPath);
    if (intp == UseBinError) // We always store variations in the bin contents
        for (int i = 0; i <= hist->GetNcells(); ++i)
            hist->SetBinContent(i, hist->GetBinError(i));
    m_global_uncs[name] = std::move(hist);
}
 
template<class Object, class... Context>
class CorrectionlibSFApplier
{
protected:
    correction::Correction::Ref m_cref; 
 
private:
    bool m_correction,    
         m_uncertainties; 
 
    std::map<std::string, correction::Variable::Type> m_bin_wise_uncs;
    std::map<std::string, correction::Variable::Type> m_global_uncs;
 
    mutable std::vector<correction::Variable::Type> m_cached_inputs;
 
public:
    CorrectionlibSFApplier (const std::filesystem::path& filePath,
                            const std::string& name,
                            bool correction,
                            bool uncertainties);
    virtual ~CorrectionlibSFApplier () noexcept = default;
 
    void operator() (Object& object, const Context&... ctx) const;
    void operator() (std::vector<Object>& objects, const Context&... ctx) const;
    std::vector<std::string> weightNames () const;
 
protected:
    void addBinWiseUnc (const std::string& name,
                        const correction::Variable::Type& input);
    void addGlobalUnc (const std::string& name,
                       const correction::Variable::Type& input);
 
    virtual bool passes (const Object& obj, const Context&... ctx) const = 0;
 
    virtual void inputs (std::vector<correction::Variable::Type>& inputs,
                         const Object& obj,
                         const Context&... ctx) const = 0;
 
    virtual std::pair<double, double>
    varied(std::vector<correction::Variable::Type>& inputs, double sf,
           const correction::Variable::Type& uncname) const = 0;
 
    virtual int binIndex (const Object& obj, const Context&... ctx) const
    { return 0; }
};
 
template<class Object, class... Context>
CorrectionlibSFApplier<Object, Context...>::CorrectionlibSFApplier (
        const std::filesystem::path& filePath,
        const std::string& name,
        bool correction,
        bool uncertainties)
    : m_correction(correction)
    , m_uncertainties(uncertainties)
{
    namespace fs = std::filesystem;
 
    if (!correction)
        return;
 
    // Check that the file exists
    if (!fs::exists(filePath))
        BOOST_THROW_EXCEPTION(
            fs::filesystem_error("Not found",
                                 filePath,
                                 make_error_code(std::errc::no_such_file_or_directory)));
 
    auto cset = correction::CorrectionSet::from_file(filePath);
    m_cref = cset->at(name);
}
 
template<class Object, class... Context>
void CorrectionlibSFApplier<Object, Context...>::operator() (
        Object& obj, const Context&... ctx) const
{
    auto& weights = obj.weights;
 
    // For safety
    if (weights.empty()) weights.emplace_back(DAS::Weight{1, 0});
 
    const auto old = weights.front().v;
 
    // Apply the cut
    if (!passes(obj, ctx...)) {
        weights *= 0;
        return; // Skip applying weights
    }
 
    // Nominal correction
    if (!m_correction) return;
 
    m_cached_inputs.clear();
    inputs(m_cached_inputs, obj, ctx...);
    const auto sf = m_cref->evaluate(m_cached_inputs);
 
    weights *= sf;
 
    // Systematics
    if (!m_uncertainties) return;
 
    // The order of the next two blocks must match weightNames()
    if (!m_bin_wise_uncs.empty()) {
        const int bin = binIndex(obj, ctx...);
        for (const auto& [_, uncname] : m_bin_wise_uncs) {
            const auto& [up, down] = varied(m_cached_inputs, sf, uncname);
            weights.push_back({old * up, bin});
            weights.push_back({old * down, bin});
        }
    }
 
    for (const auto& [_, uncname] : m_global_uncs) {
        const auto& [up, down] = varied(m_cached_inputs, sf, uncname);
        weights.push_back({old * up, 0});
        weights.push_back({old * down, 0});
    }
}
 
template<class Object, class... Context>
void CorrectionlibSFApplier<Object, Context...>::operator() (
        std::vector<Object>& objects, const Context&... ctx) const
{
    for (auto& obj: objects) {
        (*this)(obj, ctx...);
    }
}
 
template<class Object, class... Context>
std::vector<std::string>
CorrectionlibSFApplier<Object, Context...>::weightNames () const
{
    std::vector<std::string> names;
    // The order of the next two blocks must match the call operator
    for (const auto& [name, _] : m_bin_wise_uncs) {
        names.push_back(name + SysUp);
        names.push_back(name + SysDown);
    }
    for (const auto& [name, _] : m_global_uncs) {
        names.push_back(name + SysUp);
        names.push_back(name + SysDown);
    }
    return names;
}
 
template<class Object, class... Context>
void CorrectionlibSFApplier<Object, Context...>::addBinWiseUnc (
        const std::string& name, const correction::Variable::Type& input)
{
    if (!m_correction || !m_uncertainties) return;
    m_bin_wise_uncs[name] = std::move(input);
}
 
template<class Object, class... Context>
void CorrectionlibSFApplier<Object, Context...>::addGlobalUnc (
        const std::string& name, const correction::Variable::Type& input)
{
    if (!m_correction || !m_uncertainties) return;
    m_global_uncs[name] = std::move(input);
}
 
} // namespace DAS