GenericSFApplier.h

Go to the documentation of this file.

#pragma once
 
#include "Core/Objects/interface/Format.h"
#include "Core/Objects/interface/Weight.h"
 
#include <darwin.h>
 
#include <TFile.h>
#include <TH1.h>
 
#include <filesystem>
#include <map>
#include <utility>
#include <vector>
 
namespace DAS {
 
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>
Weights &weightsRef (Object& obj)
{
    return obj.weights;
}
 
template<class Object, class... Context>
class GenericSFApplier
{
    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:
    GenericSFApplier (const std::filesystem::path& filePath,
                      bool correction,
                      bool uncertainties);
    virtual ~GenericSFApplier () 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>
GenericSFApplier<Object, Context...>::GenericSFApplier (
        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 GenericSFApplier<Object, Context...>::operator() (
        Object& obj, const Context&... ctx) const
{
    auto& weights = weightsRef(obj);
 
    // For safety
    if (weights.empty()) weights.emplace_back(DAS::Weight{1, 0});
 
    const float old = weights.front().v;
    if (old == 0) return; // Don't bother
 
    // Apply the cut
    if (!passes(obj, ctx...)) {
        weights = {{0, 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 GenericSFApplier<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>
GenericSFApplier<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 GenericSFApplier<Object, Context...>::loadNominal (
        const std::string& histPath)
{
    if (!m_correction) return;
    m_nominal = getHistSafe(m_file, histPath);
}
 
template<class Object, class... Context>
void GenericSFApplier<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 GenericSFApplier<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);
}
 
} // namespace DAS