common.h

Go to the documentation of this file.

#pragma once
 
#include <cassert>
#include <iostream>
#include <optional>
#include <regex>
#include <string>
#include <sstream>
 
#include "Core/CommonTools/interface/binnings.h"
#include "Core/Objects/interface/Jet.h"
#include "Core/Objects/interface/Event.h"
 
#include <boost/algorithm/string.hpp>
#include <boost/assert/source_location.hpp>
 
#include <correction.h>
#include <darwin.h>
 
namespace DAS::JetEnergy {
 
inline static const std::vector<double> pt_JERC_edges = {15, 17, 20, 23, 27, 30, 35, 40, 45, 57, 72, 90, 120, 150, 200, 300, 400, 550, 750, 1000, 1500, 2000, 2500, 3000, 3500, 4000, 4500, 5000, 5500, 6000, 6500, 7000, 7500, 8000, 8500, 9000, 9500};
inline static const int nPtJERCbins = pt_JERC_edges.size()-1;
 
inline static const std::vector<double> abseta_edges = { 0.000, 0.261, 0.522, 0.783, 1.044, 1.305, 1.566, 1.740, 1.930, 2.043, 2.172, 2.322, 2.500, 2.650, 2.853, 2.964, 3.139, 3.489, 3.839, 5.191 };
inline static const int nAbsEtaBins = abseta_edges.size()-1;
inline static const std::vector<TString> absetaBins = MakeTitle(abseta_edges, "|#eta|", false, true, [](double v) { return Form("%.3f", v);});
 
inline static const std::map<int,std::vector<double>> rho_edges = {
    { 2016, {0, 6.69, 12.39, 18.09, 23.79, 29.49, 35.19, 40.9, 999} }, // EOY16
    { 2017, {0, 7.47, 13.49, 19.52, 25.54, 31.57, 37.59, 999} }, // UL17
    { 2018, {0, 7.35, 13.26, 19.17, 25.08, 30.99, 36.9, 999} }, // UL18
 
    { 2024, {0.0, 8.23, 14.49, 20.75, 27.01, 33.27, 39.53, 45.79, 52.05, 999} }, // Winter24 (copy of Summer23)
};
 
inline int GetNRhoBins (int year) { return rho_edges.at(year).size()-1; }
inline std::vector<TString> rhoBins (int year)
{
    return MakeTitle(rho_edges.at(year), "#rho", false, false,
                        [](double v) { return Form("%.2f", v); });
}
 
inline static const std::vector<double> eta_unc_edges = {-5.4, -5.0, -4.4, -4.0, -3.5, -3.0, -2.8, -2.6, -2.4, -2.2, -2.0, -1.8, -1.6, -1.4, -1.2, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0 , 0.2 , 0.4 , 0.6 , 0.8 , 1.0 , 1.2 , 1.4 , 1.6 , 1.8 , 2.0 , 2.2 , 2.4 , 2.6 , 2.8 , 3.0, 3.5 , 4.0 , 4.4 , 5.0 , 5.4};
inline static const int nEtaUncBins = eta_unc_edges.size()-1;
 
inline static const std::vector<std::string> JES_variations {"AbsoluteStat", "AbsoluteScale",  "AbsoluteMPFBias", "Fragmentation", "SinglePionECAL", "SinglePionHCAL", "FlavorQCD", "RelativeJEREC1", "RelativeJEREC2", "RelativeJERHF", "RelativePtBB", "RelativePtEC1", "RelativePtEC2", "RelativePtHF", "RelativeBal", "RelativeSample", "RelativeFSR", "RelativeStatFSR", "RelativeStatEC", "RelativeStatHF", "PileUpDataMC", "PileUpPtRef", "PileUpPtBB", "PileUpPtEC1", "PileUpPtEC2", "PileUpPtHF"}; 
 
inline static const float w = 0.8;
 
inline std::vector<double> getBinning (int nBins, float first, float last)
{
    std::vector<double> bins(nBins+1);
    for(int i = 0; i <= nBins; ++i)
        bins[i] = first + ((last-first)/nBins) * i;
    return bins;
}
 
inline void assertValidBinning (const std::vector<double>& v)
{
    assert(v.size() > 1);
    for (size_t i = 1; i<v.size(); ++i) {
        if (v[i] > v[i-1]) continue;
        std::cerr << i << ' ' << v[i] << ' ' << v[i-1] << '\n';
    }
}
 
inline std::vector<float> GetLogBinning (float minpt, float maxpt, int nbins)
{
    assert(maxpt > minpt);
    assert(minpt > 0);
    assert(nbins > 1);
 
    std::vector<float> edges;
 
    float R = pow(maxpt/minpt,1./nbins);
    for (float pt = minpt; pt <= maxpt; pt *= R) edges.push_back(pt);
    std::cout << edges.size() << std::endl;
 
    return edges;
}
 
inline std::string GetAlgo (const Darwin::Tools::UserInfo& metainfo)
{
    if (metainfo.Find("flags", "labels", "CHS"  )) return "chs";
    if (metainfo.Find("flags", "labels", "PUPPI")) return "Puppi";
 
    std::cerr << orange << "Couldn't identify CHS or PUPPI. Running default "
                           "CHS for Run 2 (PUPPI foir Run 3)\n" << def;
 
    auto year = metainfo.Get<int>("flags", "year");
    if (year > 2019) return "Puppi"; // run 3
    if (year > 2014) return "chs"; // run 2
                     return ""; // run 1
}
 
std::string GetShortCampaign (const std::string& campaign)
{
    using namespace std;
    regex r("^(Summer|Fall|Autumn|Winter|Spring)[0-9]{2}[A-Za-z0-9]*");
    smatch campaign_short;
    if (!regex_search(campaign, campaign_short, r))
        BOOST_THROW_EXCEPTION( invalid_argument("The campaign could not be identified"));
    return campaign_short.str();
}
 
std::string ScanCorrections (const auto& corrections, 
                             const std::string& key = "") 
{
    std::stringstream s;
    s << "Available corrections:";
    for (const auto& correction: corrections) {
        bool found = correction.first == key;
        if (found) s << bold << green;
        s << ' ' << correction.first;
        if (found) s << def;
    }
    return s.str();
}
 
template<typename CorrectionType> 
CorrectionType GetCorrection (const auto& corrections, 
                              const std::string& campaign, 
                              const std::string& type, 
                              const std::string& level, 
                              const std::string& suffix) 
{
    namespace al = boost::algorithm;
    for (const auto& correction: corrections) {
        if (!al::starts_with(correction.first, campaign) ||
            !al::contains   (correction.first, type) ||
            !al::contains   (correction.first, level) ||
            !al::ends_with  (correction.first, suffix))
            continue;
 
        return correction.second;
    }
    using namespace std;
    BOOST_THROW_EXCEPTION(
            invalid_argument(
                "No `"s + campaign + "*"s + type + "*" + level + "*"s + suffix
                + "` correction can be found in the given tables. "s
                + ScanCorrections(corrections)
            )
        );
}
 
float Evaluate (const auto& correction, 
                std::ostream& cout, 
                const RecJet& recJet, 
                const std::optional<float>& rho = {}, 
                const std::optional<GenJet>& genJet = {}, 
                const std::optional<RecEvent>& recEvt = {}, 
                const std::optional<std::string>& systematic = {}, 
                const std::optional<float>& JER = {}, 
                const std::optional<float>& JERSF = {}  
               )
try {
    using namespace std;
 
    cout << correction->name();
    vector<correction::Variable::Type> inputs;
    inputs.reserve(10);
 
    auto push_back = [&cout,&inputs](const auto& v) {
        inputs.push_back(v);
        cout << '=' << v;
    };
 
    for (const correction::Variable& input: correction->inputs()) {
        const string& n = input.name();
        cout << ' ' << n;
        if      (n == "JetPt"     ) push_back(recJet.CorrPt());
        else if (n == "JetEta"    ) push_back(recJet.p4.Eta());
        else if (n == "JetPhi"    ) push_back(recJet.p4.Phi());
        else if (n == "JetEta"    ) push_back(recJet.p4.Phi());
        else if (n == "JetA"      ) push_back(recJet.area);
        else if (n == "Rho"       ) push_back(*rho);
        else if (n == "systematic") push_back(*systematic);
        else if (n == "GenPt"     ) push_back(genJet ? genJet->CorrPt() : -1);
        else if (n == "EventID"   ) push_back((int) recEvt->evtNo);
        else if (n == "run"       ) push_back((float) recEvt->runNo);
        else if (n == "JER"       ) push_back(*JER);
        else if (n == "JERSF"     ) push_back(*JERSF);
        else BOOST_THROW_EXCEPTION( invalid_argument("`"s + n + "` is needed by "s
                                    + correction->name() + " but not recognized."s) );
    }
    auto corr = correction->evaluate(inputs);
    cout << " -> " << corr << endl;
    return corr;
}
catch (std::runtime_error& e) {
    const char * name = typeid(*correction).name();
    auto location = boost::source_location(__FILE__, __LINE__, name);
    boost::throw_exception(e, location);
}
 
int GetR (const Darwin::Tools::UserInfo& metainfo)
{
    auto year = metainfo.Get<int>("flags", "year"),
         R = metainfo.Get<int>("flags", "R");
 
    std::string warning = "Not a standard jet size.\n";
    if (year > 2014) { // run 2 and run 3
        if (R != 4 && R != 8)
            std::cerr << orange << warning << def;
        return R < 6 ? 4 : 8;
    }
    else { // run 1
        if (R != 5 && R != 7)
            std::cerr << orange << warning << def;
        return R < 6 ? 5 : 7;
    }
}
 
} // end of DAS::JetEnergy namespace