Balance.h

Go to the documentation of this file.

#pragma once
 
#include <map>
#include <memory>
#include <vector>
 
#include <TDirectory.h>
#include <TH3.h>
#include <TList.h>
#include <TString.h>
 
#include <darwin.h>
 
#include "Core/Objects/interface/Jet.h"
#include "Core/Objects/interface/Variation.h"
 
#include "Core/JEC/interface/common.h"
#include "Core/CommonTools/interface/toolbox.h"
 
namespace DAS::JetEnergy {
 
using namespace std;
namespace DT = Darwin::Tools;
 
vector<DAS::Uncertainties::Variation> GetScaleVariations (const DT::MetaInfo& metainfo,
                                                          bool applySyst)
{
    vector<DAS::Uncertainties::Variation> variations;
    variations.emplace_back("nominal","nominal");
    if (!applySyst) return variations;
 
    const TList * groupContents = metainfo.List("variations", DAS::RecJet::ScaleVar);
    if (!groupContents)
        BOOST_THROW_EXCEPTION( invalid_argument(groupContents->GetName()) );
 
    size_t i = 0;
    for (const TObject * obj: *groupContents) {
        const auto os = dynamic_cast<const TObjString*>(obj);
        if (!os) BOOST_THROW_EXCEPTION( invalid_argument(obj->GetName()) );
        TString s = os->GetString();
        s.ReplaceAll(SysUp, "");
        s.ReplaceAll(SysDown, "");
        if (find(begin(JES_variations), end(JES_variations), s) != end(JES_variations))
            variations.emplace_back(DAS::RecJet::ScaleVar, os->GetString(), ++i);
    }
    return variations;
}
 
class Balance {
 
    TDirectory * d; 
    unique_ptr<TH3D> asymmetry; 
    unique_ptr<TH2D> response; 
 
public:
    DAS::Uncertainties::Variation v; 
 
    Balance (TDirectory * d,
             const DAS::Uncertainties::Variation& v,
             const vector<double>& pt_edges,
             const vector<double>& y_edges) :
        d(d), v(v)
    {
        int nBins = 120;
        vector balBins = MakeBinning<double>({-1.2,1.2}, nBins);
 
        int nPtBins = pt_edges.size()-1,
            nYbins  =  y_edges.size()-1;
 
        TString dirname = d->GetName();
        TString hname = dirname + '_' + v.group + '_' + v.name;
        asymmetry = make_unique<TH3D>(hname, v.name, nBins, balBins.data(),
                                                     nPtBins, pt_edges.data(),
                                                     nYbins, y_edges.data());
 
        balBins = MakeBinning<double>({0., 1.2}, nBins);
        response = make_unique<TH2D>("response_" + hname, v.name, nBins, balBins.data(),
                                                                  nYbins, y_edges.data());
    }
 
    void operator() (const DAS::AbstractPhysicsObject& tag,
                     const DAS::AbstractPhysicsObject& probe,
                     const DAS::MET* met,
                     const double ev_w)
    {
        float absy_probe = probe.AbsRap(v);
 
        float pt_tag   = tag  .CorrPt(v),
              pt_probe = probe.CorrPt(v);
 
        double jw = tag  .Weight(v)
                  * probe.Weight(v);
 
        response->Fill(pt_probe/pt_tag, absy_probe, ev_w*jw);
 
        // \todo does this also apply to Z+jet balance?
        //       (e.g. not in 1804.05252)
        float pt_dif = (pt_tag - pt_probe),
              pt_ave = (pt_tag + pt_probe) * 0.5;
 
        float asym = pt_dif / pt_ave; // pT balance method
 
        if (met != nullptr) // MPF method ref:1607.03663 $\mathcal{A} =\frac{\vec{p}_T^{miss}\cdot\left(\vec{p}_{T,tag}/p_{T,tag}\right)}{2 p_{T,ave}}$
            asym = -met->Pt * cos(met->Phi - tag.CorrP4(v).Phi()) / pt_ave;
 
        asymmetry->Fill(asym, pt_ave, absy_probe, ev_w*jw);
    }
 
    void Write ()
    {
        d->cd();
        TString dirname = d->GetName();
 
        asymmetry->SetDirectory(d);
        TString hname = asymmetry->GetName();
        hname.ReplaceAll(dirname + '_', "");
        asymmetry->Write(hname);
 
        response->SetDirectory(d);
        hname = response->GetName();
        hname.ReplaceAll(dirname + '_', "");
        response->Write(hname);
    }
};
 
} // end of DAS::JetEnergy namespace