Plots.h

Go to the documentation of this file.

#pragma once
 
#include <algorithm>
#include <memory>
#include <vector>
 
#include <TH2F.h>
 
#include "Core/Objects/interface/Di.h"
#include "Core/Objects/interface/Event.h"
#include "Core/Objects/interface/Jet.h"
 
#include "Core/CommonTools/interface/binnings.h"
 
namespace DAS::PUstaub {
 
static const std::vector<double> pthat_edges = {0, 1e-6, 15, 30, 50, 80, 120, 170, 300, 470, 600, 800, 1000, 1400, 1800, 2400, 3200, 6500};
static const int nPtHatBins = pthat_edges.size()-1;
 
struct Plots {
 
    TString name;
 
    std::unique_ptr<TH1F> pthatME, pthatPU, pthatMax, genht, recht;
    std::unique_ptr<TH2F> pthatLogw, genptLogw, genhtLogw, genmjjLogw, recptLogw, rechtLogw, recmjjLogw;
 
    Plots (const char * Name) :
        name(Name),
        pthatME(std::make_unique<TH1F>("pthatME", ";#hat{p}_{T}^{ME}   [GeV];N_{eff}^{collisions}", nPtHatBins, pthat_edges.data())),
        pthatPU(std::make_unique<TH1F>("pthatPU", ";#hat{p}_{T}^{PU}   [GeV];N_{eff}^{collisions}", nPtHatBins, pthat_edges.data())),
        pthatMax(std::make_unique<TH1F>("pthatMax", ";max #hat{p}_{T}^{ME+PU}   [GeV];N_{eff}^{bx}", nPtHatBins, pthat_edges.data())),
        genht(std::make_unique<TH1F>("genht", ";H_{T}^{gen}   [GeV];N_{eff}^{j}", nPtBins, pt_edges.data())),
        recht(std::make_unique<TH1F>("recht", ";H_{T}^{rec}   [GeV];N_{eff}^{j}", nPtBins, pt_edges.data())),
        pthatLogw(std::make_unique<TH2F>("pthatLogw", ";#hat{p}_{T}   [GeV];log(w);N_{eff}^{j}", nPtHatBins, pthat_edges.data(), 400, -30, 20)),
        genptLogw(std::make_unique<TH2F>("genptLogw", ";Jet p_{T}^{gen}   [GeV];log(w);N_{eff}^{j}", nPtBins, pt_edges.data(), 400, -30, 20)),
        genhtLogw(std::make_unique<TH2F>("genhtLogw", ";H_{T}^{gen}   [GeV];log(w);N_{eff}^{j}", nPtBins, pt_edges.data(), 400, -30, 20)),
        genmjjLogw(std::make_unique<TH2F>("genmjjLogw", ";m_{jj}^{gen}   [GeV];log(w);N_{eff}^{jj}", nMjjBins, Mjj_edges.data(), 400, -30, 20)),
        recptLogw(std::make_unique<TH2F>("recptLogw", ";Jet p_{T}^{rec}   [GeV];log(w);N_{eff}^{j}", nPtBins, pt_edges.data(), 400, -30, 20)),
        rechtLogw(std::make_unique<TH2F>("rechtLogw", ";H_{T}^{rec}   [GeV];log(w);N_{eff}^{j}", nPtBins, pt_edges.data(), 400, -30, 20)),
        recmjjLogw(std::make_unique<TH2F>("recmjjLogw", ";m_{jj}^{rec}   [GeV];log(w);N_{eff}^{jj}", nMjjBins, Mjj_edges.data(), 400, -30, 20))
    { }
 
private:
 
    float FillPtHat (const GenEvent& event, const PileUp& pileup)
    {
        float gW = event.weights.front();
 
        pthatME->Fill(event.hard_scale, gW);
        float maxpthat = event.hard_scale;
        for (auto pthat: pileup.pthats) {
            pthatPU->Fill(pthat, gW);
            maxpthat = std::max(maxpthat, pthat);
        }
        pthatMax->Fill(maxpthat, gW);
 
        return maxpthat;
    }
 
    template<typename Jet> void FillObsWgt (std::unique_ptr<TH2F>& ptLogw,
                                            std::unique_ptr<TH2F>& mjjLogw,
                                            std::unique_ptr<TH2F>& htLogw,
                                            std::unique_ptr<TH1F>& ht,
                                            const std::vector<Jet>& jets, float evW)
    {
        auto pt_sum = 0;
        auto weight_sum = 0;
        for (const auto& jet: jets) {
            auto jW = log(jet.weights.front());
            weight_sum = weight_sum + jW;
            ptLogw->Fill(jet.p4.Pt(), evW + jW);
            pt_sum += jet.p4.Pt();
        }
        ht->Fill(pt_sum);
 
        htLogw->Fill(pt_sum, evW + weight_sum);
 
        if (jets.size() < 2) return;
        const auto& j0 = jets.at(0),
                    j1 = jets.at(1);                                     
 
        if (j0.p4.Pt() < 100. || j1.p4.Pt() < 50.) return;
        if (std::max(j0.AbsRap(), j1.AbsRap()) >= 3.0) return;
        auto dijet = j0 + j1;
        auto djW = log(j0.weights.front()) + log(j1.weights.front());
        mjjLogw->Fill(dijet.CorrP4().M(), evW + djW);
    }
 
public:
 
    void operator() (const std::vector<GenJet>& genjets,
                     const std::vector<RecJet>& recjets,
                     const GenEvent& genEvt,
                     const RecEvent& recEvt,
                     const PileUp& pileup)
    {
        float maxpthat = FillPtHat(genEvt, pileup);
 
        float logw = log(genEvt.weights.front());
        pthatLogw->Fill(maxpthat, logw);
        FillObsWgt<GenJet>(genptLogw, genmjjLogw, genhtLogw, genht, genjets, logw);
 
        logw += log(recEvt.weights.front());
        FillObsWgt<RecJet>(recptLogw, recmjjLogw, rechtLogw, recht, recjets, logw);
 
    }
 
    void Write (TDirectory * d)
    {
        d->cd();
        TDirectory * dd = d->mkdir(name);
        dd->cd();
        std::vector<TH1*> hists {pthatME.get(), pthatPU.get(), pthatMax.get(), genht.get(), recht.get(), 
                            dynamic_cast<TH1*>(pthatLogw.get()),
                            dynamic_cast<TH1*>(genptLogw.get()),
                            dynamic_cast<TH1*>(genhtLogw.get()),
                            dynamic_cast<TH1*>(genmjjLogw.get()),
                            dynamic_cast<TH1*>(recptLogw.get()),
                            dynamic_cast<TH1*>(rechtLogw.get()),
                            dynamic_cast<TH1*>(recmjjLogw.get())};
        for (auto h: hists) {
            h->SetDirectory(dd);
            TString hname = h->GetName();
            hname.ReplaceAll(name, "");
            h->Write(hname);
        }
    }
};
 
}