TriggerEff.h

Go to the documentation of this file.

#pragma once
 
#include <cstdlib>
 
#include <TH2.h>
#include <TFile.h>
 
#include <filesystem>
#include <iostream>
#include <map>
#include <memory>
#include <vector>
 
#include "Core/CommonTools/interface/binnings.h"
#include "Core/Normalisation/interface/TriggerLumi.h"
#include "Core/Objects/interface/Jet.h"
#include "Core/Trigger/interface/match.h"
 
#include <darwin.h>
 
namespace fs = std::filesystem;
 
namespace DAS::Normalisation {
 
struct TriggerEff {
 
    TH2 * h; 
 
    std::vector<TH2 *> contribs; 
 
    TriggerEff (const fs::path& file, const std::map<int, TriggerLumi>& triggers_lumi, int year) :
        h(new TH2D("efficiency", "efficiency", nPtBins, pt_edges.data(), nYbins, y_edges.data()))
    {
        using namespace std;
 
        assert(fs::exists(file)); 
        auto f = make_unique<TFile>(file.c_str(), "READ");
 
        int firsttrigger = triggers_lumi.begin()->first; // e.g. typically 40 for HLT_PFJet_* in high pile-up runs
 
        for (auto trigger_lumi: triggers_lumi) {
            int trigger = trigger_lumi.first,
                turnon = trigger_lumi.second.turnon;
 
            const char * method = trigger == firsttrigger ? "TnP" : "emulation";
            auto hIn = unique_ptr<TH2>(f->Get<TH2>(Form("%s/HLT%d/efficiency", method, trigger)));
            if (hIn.get() == nullptr) {
                cerr << red << "No efficiency found for trigger " << trigger << normal << '\n';
                contribs.push_back(nullptr);
                continue;
            }
            cout << "Getting efficiency for trigger " << trigger << ": " << hIn.get() << endl;
 
            assert(hIn->GetNbinsX() == h->GetNbinsX());
 
            int I = h->GetXaxis()->FindBin(turnon);
 
            for (int iy = 1; iy <= nYbins; ++iy)
            for (int ipt = I; ipt <= h->GetNbinsX(); ++ipt) {
                double content = hIn->GetBinContent(ipt, iy),
                       error   = hIn->GetBinError  (ipt, iy);
                h->SetBinContent(ipt, iy, content);
                h->SetBinError  (ipt, iy, error  );
            }
 
            auto name = Form("HLT%d", trigger);
            auto hContrib = dynamic_cast<TH2*>(hIn->Clone(name));
            hContrib->SetDirectory(0);
            hContrib->Reset();
            contribs.push_back(hContrib);
        }
        f->Close();
 
        // regularise at high pt
        auto lasttrigger = prev(triggers_lumi.end())->first;
        for (int iy = 1; iy <= nYbins; ++iy) {
            double turnon = triggers_lumi.at(lasttrigger).turnon;
            int ipt = h->GetXaxis()->FindBin(turnon);
 
            // first, find the first bin of the last trigger reaching 0.9995
            while (ipt <= h->GetNbinsX()) {
                double content = h->GetBinContent(ipt, iy);
                if (content > 0.9995) break; 
                ++ipt;
            }
            // then, put all next bins to one
            while (ipt <= h->GetNbinsX()) {
                h->SetBinContent(ipt,iy,1);
                h->SetBinError  (ipt,iy,0);
                ++ipt;
            }
        }
    }
 
    double operator()(const DAS::RecJet& j)
    {
        using namespace std;
 
        auto pt = j.CorrPt();
        if (pt > 905) return 1.; 
 
        auto y = j.AbsRap();
 
        int ipt = h->GetXaxis()->FindBin(pt),
            iy  = h->GetYaxis()->FindBin(y );
 
        auto e = h->GetBinContent(ipt, iy);
        if (e < 0 || e > 1) {
            cerr << "\x1B[31m\e[1m" << pt << ' ' << y << '\t' << ipt << ' ' << iy << '\t' << e << '\n';
            h->Print("all");
            cerr << "\x1B[30m\e[0m";
            exit(EXIT_FAILURE);
        }
        return e;
    }
};
 
} // end of DAS::Normalisation namespace