PUcorrection.h

Go to the documentation of this file.

#include <cstdlib>
#include <cassert>
#include <iostream>
#include <limits>
#include <functional>
 
#include <TFile.h>
#include <TH1.h>
 
#include <Core/Objects/interface/Format.h>
 
#include <darwin.h>
 
namespace DAS::PUprofile {
 
struct Correction {
    TH1 * nominal, * upper, * lower;
 
    Correction (TFile * fMC, TFile * fData,
                std::function<const char *(const char *, bool)> hname,
                float maxWeight = std::numeric_limits<float>::max()
               ) :
        nominal(fData->Get<TH1>(hname("nominal", false))),
        upper  (fData->Get<TH1>(hname(("PU" + SysUp).c_str(), false))),
        lower  (fData->Get<TH1>(hname(("PU" + SysDown).c_str(), false)))
    {
        namespace DE = Darwin::Exceptions;
        for (TH1 * h: {nominal, upper, lower}) {
            if (h == nullptr)
                BOOST_THROW_EXCEPTION( DE::BadInput("No profile found in input", *fData) );
            h->SetDirectory(0);
            auto integral = h->Integral();
            if (integral == 0)
                BOOST_THROW_EXCEPTION( DE::BadInput("Empty profile", *h) );
            h->Scale(1./integral);
        }
 
        auto n = hname("nominal", true);
        std::cout << "Fetching " << n << std::endl;
        auto sim = dynamic_cast<TH1*>(fMC->Get(n));
        if (sim == nullptr)
            BOOST_THROW_EXCEPTION( DE::BadInput("No profile found in input", *fMC) );
        auto integral = sim->Integral();
        if (integral == 0)
            BOOST_THROW_EXCEPTION( DE::BadInput("Empty profile", *sim) );
        sim->Scale(1./integral);
        for (TH1 * h: {nominal, upper, lower})
            h->Divide(sim);
 
        for (TH1 *h: {nominal, upper, lower}) {
            for (int i = 0; i <= h->GetNbinsX()+1; ++i) {
                if (h->GetBinContent(i) < maxWeight) continue;
                h->SetBinContent(i, maxWeight);
            }
        }
    }
 
    float operator() (float pu, char v = '0') const
    {
        namespace DE = Darwin::Exceptions;
        TH1 * h = nullptr;
        switch (v) {
            case '0': h = nominal; break;
            case '+': h = upper  ; break;
            case '-': h = lower  ; break;
            default:
                      BOOST_THROW_EXCEPTION( std::runtime_error("Unknown variation") );
        }
        if (h == nullptr)
            BOOST_THROW_EXCEPTION( DE::BadInput("No correction found for this variation input", *h) );
        int ibin = h->FindBin(pu);
        return h->GetBinContent(ibin);
    }
 
    void Write (TDirectory * d, bool reset)
    {
        d->cd();
        for (TH1 * h: {nominal, upper, lower}) {
            h->SetDirectory(d);
            if (reset) h->Reset();
        }
        nominal->Write("nominal");
        upper->Write("upper");
        lower->Write("lower");
    }
};
 
} // end of namespaces