applyPrefiringWeights.h

Go to the documentation of this file.

#pragma once
 
#include <cassert>
#include <cstdlib>
 
#include <iostream>
#include <filesystem>
#include <tuple>
#include <vector>
 
#include "Core/Objects/interface/Event.h"
#include "Core/Objects/interface/Jet.h"
 
#include "Core/CommonTools/interface/binnings.h"
#include "Core/CommonTools/interface/toolbox.h"
 
#include <TF1.h>
#include <TH2.h>
 
#include <darwin.h>
 
namespace fs = std::filesystem;
namespace DE = Darwin::Exceptions;
 
static const auto feps = std::numeric_limits<float>::epsilon();
 
namespace DAS::Prefiring {
 
struct Functor {
 
    std::tuple<TH2*, TH2*, TH2*> prefBinnedWeights; 
    std::tuple<TF1 *, TF1*> prefSmoothWeights; 
 
    using Operation = std::tuple<float,float,float> (Functor::*)(int,const RecJet&) const;
    const Operation GetWeights;
 
    const bool syst;
 
    std::vector<TH2 *> multiplicities; 
 
    static std::tuple<TH2*, TH2*, TH2*> PrepareBinnedWeights (const std::unique_ptr<TH2>& hIn)
    {
        using namespace std;
 
        auto MapNo = dynamic_cast<TH2*>(hIn->Clone("nominal_map")),
             MapUp = dynamic_cast<TH2*>(hIn->Clone("upper_map"  )),
             MapDn = dynamic_cast<TH2*>(hIn->Clone("lower_map"  ));
        MapNo->SetDirectory(0);
        MapUp->SetDirectory(0);
        MapDn->SetDirectory(0);
        for (int etabin = 0; etabin <= hIn->GetNbinsX()+1; ++etabin)
        for (int  ptbin = 0;  ptbin <= hIn->GetNbinsY()+1; ++ ptbin) {
            auto proba = hIn->GetBinContent(etabin, ptbin),
                 error = hIn->GetBinError  (etabin, ptbin);
 
            // uncertainty only in the forward region! (the following boolean is a trick :D)
            bool Forward = (proba > feps); // means proba > 0
            auto unc = Forward ? max(proba * 0.2f, error) : 0.f;
 
            // the map contains the probability, but we want the weight
            auto weight = 1. - proba;
 
            MapNo->SetBinContent(etabin, ptbin, weight);
            MapUp->SetBinContent(etabin, ptbin, weight+unc);
            MapDn->SetBinContent(etabin, ptbin, weight-unc);
 
            //if (Forward) cout << green;
            //cout << setw(10) << etabin << setw(10) << ptbin << setw(20) << proba << setw(15) << error << setw(20) << weight << setw(15) << unc << /*normal <<*/ endl;
        }
 
        return { MapNo, MapUp, MapDn };
    }
 
    static std::tuple<TF1*, TF1*> PrepareSmoothWeights (TString hname, std::unique_ptr<TFile> f)
    {
        using namespace std;
 
        cout << hname << endl;
        auto f5 = f->Get<TF1>(hname + "_ybin5"),
             f6 = f->Get<TF1>(hname + "_ybin6");
        assert(f5 != nullptr); 
        assert(f6 != nullptr); 
        return {f5, f6};
    }
 
    static std::vector<TH2 *> GetMultiplicityPlots (TH2 * templ, int n = 7)
    {
        using namespace std;
 
        vector<TH2*> v;
        for (int i = 1; i <= n; ++i) {
            auto h = dynamic_cast<TH2*>(templ->Clone(Form("multiplicity_nbin%d", i)));
            h->Reset();
            v.push_back(h);
        }
        return v;
    }
 
    Functor
        (const fs::path& file, 
         const std::string name, 
         const std::string type, 
         bool applySyst
         ) :
        GetWeights(type == "smooth" ? &Functor::GetSmoothWeights :
                   type == "binned" ? &Functor::GetBinnedWeights : nullptr),
        syst(applySyst)
    {
        using namespace std;
 
        if (!fs::exists(file))
            BOOST_THROW_EXCEPTION( fs::filesystem_error("File does not exists.",
                                   file, make_error_code(errc::no_such_file_or_directory)) );
 
        auto f = make_unique<TFile>(file.c_str(), "READ");
        string hname = file.stem().string();
        auto h = unique_ptr<TH2>(f->Get<TH2>(hname.c_str()));
        if (h.get() == nullptr) {
            string what = "Map " + hname + " can't be found";
            BOOST_THROW_EXCEPTION( DE::BadInput(what.c_str(), f) );
        }
        // this map contains the probabilities, not the weights (p + w = 1)
 
        // retrieve map in histograms
        if (type == "binned")
            prefBinnedWeights = PrepareBinnedWeights(h);
        else if (type == "smooth")
            prefSmoothWeights = PrepareSmoothWeights(file.c_str(), move(f));
        else
            BOOST_THROW_EXCEPTION( invalid_argument("Types may only be 'smooth' and 'binned'.") );
 
        TH2 * templ = get<0>(prefBinnedWeights);
        multiplicities = GetMultiplicityPlots(templ);
 
        cout << "Prefiring weights are ready" << endl;
    }
 
private:
    std::tuple<float,float,float> GetBinnedWeights (int run, const RecJet& recjet) const
    {
        using namespace std;
 
        TH2 * MapNo = get<0>(prefBinnedWeights),
            * MapUp = get<1>(prefBinnedWeights),
            * MapDn = get<2>(prefBinnedWeights);
 
        FourVector p4 = recjet.p4;
        p4 *= recjet.scales.front();
 
        auto eta = p4.Eta();
 
        if (std::abs(eta) < 2.0 || std::abs(eta) >= 3.0) return { 1, 1, 1};
 
        auto pt = p4.Pt();
 
        int etabin =     MapNo->GetXaxis()->FindBin(eta),
            ptbin  = min(MapNo->GetYaxis()->FindBin(pt ), MapNo->GetNbinsY()); // i.e. don't take the overflow
 
        auto corr  = MapNo->GetBinContent(etabin, ptbin),
             varUp = MapUp->GetBinContent(etabin, ptbin),
             varDn = MapDn->GetBinContent(etabin, ptbin);
 
        //cout << run << setw(15) << eta << setw(15) << pt << setw(10) << etabin << setw(10) << ptbin << setw(20) << corr << setw(10) << varUp << setw(10) << varDn << '\n';
 
        return {corr, varUp, varDn};
    }
 
    std::tuple<float, float, float> GetSmoothWeights (int run, const RecJet& recjet) const
    {
        using namespace std;
 
        FourVector p4 = recjet.p4;
        p4 *= recjet.scales.front();
 
        auto eta = p4.Eta();
        int etabin = 1 + abs(eta)*2;
 
        //cout << eta << ' ' << etabin << endl;
        if (etabin < 5 || etabin > 6) return { 1, 1, 1};
 
        auto f = etabin == 5 ? get<0>(prefSmoothWeights)
                             : get<1>(prefSmoothWeights);
        auto pt = p4.Pt();
        auto eval = f->Eval(pt);
 
        float corr  = 1-    eval,
              varUp = 1-1.5*eval,
              varDn = 1-0.75*eval;
 
        //cout << run << setw(15) << eta << setw(15) << pt << setw(10) << etabin << setw(20) << corr << setw(10) << varUp << setw(10) << varDn << '\n';
 
        return {corr, varUp, varDn};
    }
 
public:
    void operator() (RecEvent & evnt, const std::vector<RecJet>& recjets)
    {
        using namespace std;
        //cout << "=====" << endl;
 
        // first, get the weights as a function of all jets present at |eta| > 2.0
        tuple<float,float,float> wAllJets { 1, 1, 1 };
        for (const RecJet& recjet: recjets) {
 
            auto wJet = (this->*GetWeights)(evnt.runNo, recjet);
 
            get<0>(wAllJets) *= get<0>(wJet);
            get<1>(wAllJets) *= get<1>(wJet);
            get<2>(wAllJets) *= get<2>(wJet);
        }
 
        // then fill control plots
        size_t i = 0;
        for (const RecJet& recjet: recjets) {
 
            float pt = recjet.CorrPt(0),
                  wEvt0 = evnt.weights.front(),
                  wJet0 = recjet.weights.front();
 
            // multiplicites
            static const size_t n = multiplicities.size();
            if (i >= n) continue;
            multiplicities.at(i)->Fill(recjet.p4.Eta(), pt, wEvt0 * wJet0);
            ++i;
        }
 
        //cout << evnt.runNo << setw(20) << get<0>(wAllJets) << setw(20) << get<1>(wAllJets) << setw(20) << get<2>(wAllJets) << '\n';
        //cout << "-----" << endl;
 
        // finally, modify the branches
        float rw0 = evnt.weights.front();
        assert(get<0>(wAllJets) > 0); // this could only arrive for jets < 30 GeV
        evnt.weights /= get<0>(wAllJets);
 
        if (!syst) return;
        evnt.weights.push_back(Weight{rw0 / get<1>(wAllJets)});
        evnt.weights.push_back(Weight{rw0 / get<2>(wAllJets)});
    }
 
    void Write (TDirectory * dir)
    {
        using namespace std;
 
        dir->cd();
        TH2 * MapNo = get<0>(prefBinnedWeights),
            * MapUp = get<1>(prefBinnedWeights),
            * MapDn = get<2>(prefBinnedWeights);
        for (TH2 * h: {MapNo, MapUp, MapDn}) {
            if (h == nullptr) continue;
            h->SetDirectory(dir);
            TString name = h->GetName();
            name = name(0, name.First('_')-1);
            h->Write(name);
        }
        TF1 * f5 = get<0>(prefSmoothWeights),
            * f6 = get<1>(prefSmoothWeights);
        if (f5 != nullptr) f5->Write("ybin5");
        if (f6 != nullptr) f6->Write("ybin6");
 
        for (auto m: multiplicities) {
            m->SetDirectory(dir);
            m->Write();
        }
    }
};
 
} // end of namespace