FFTJetPileupEstimator.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    RecoJets/FFTJetProducers
// Class:      FFTJetPileupEstimator
// 
//
// Original Author:  Igor Volobouev
//         Created:  Wed Apr 20 13:52:23 CDT 2011
// $Id: FFTJetPileupEstimator.cc,v 1.3 2012/11/21 03:13:26 igv Exp $
//
//

#include <cmath>

// Framework include files
#include "FWCore/Framework/interface/EDProducer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"

// Data formats
#include "DataFormats/Common/interface/View.h"
#include "DataFormats/Common/interface/Handle.h"
// #include "DataFormats/Histograms/interface/MEtoEDMFormat.h"
#include "DataFormats/JetReco/interface/FFTJetPileupSummary.h"
#include "DataFormats/JetReco/interface/DiscretizedEnergyFlow.h"

#include "RecoJets/FFTJetProducers/interface/FFTJetParameterParser.h"

// Loader for the lookup tables
#include "JetMETCorrections/FFTJetObjects/interface/FFTJetLookupTableSequenceLoader.h"

#define init_param(type, varname) varname (ps.getParameter< type >( #varname ))

using namespace fftjetcms;

//
// class declaration
//
class FFTJetPileupEstimator : public edm::EDProducer
{
public:
    explicit FFTJetPileupEstimator(const edm::ParameterSet&);
    ~FFTJetPileupEstimator();

protected:
    // methods
    void beginJob();
    void produce(edm::Event&, const edm::EventSetup&);
    void endJob();

private:
    FFTJetPileupEstimator();
    FFTJetPileupEstimator(const FFTJetPileupEstimator&);
    FFTJetPileupEstimator& operator=(const FFTJetPileupEstimator&);

    std::auto_ptr<reco::FFTJetPileupSummary> calibrateFromConfig(
        double uncalibrated) const;

    std::auto_ptr<reco::FFTJetPileupSummary> calibrateFromDB(
        double uncalibrated, const edm::EventSetup& iSetup) const;

    template<class Ptr>
    inline void checkConfig(const Ptr& ptr, const char* message)
    {
        if (ptr.get() == NULL)
            throw cms::Exception("FFTJetBadConfig") << message << std::endl;
    }

    edm::InputTag inputLabel;
    std::string outputLabel;
    double cdfvalue;
    double ptToDensityFactor;
    unsigned filterNumber;
    std::vector<double> uncertaintyZones;
    std::auto_ptr<fftjet::Functor1<double,double> > calibrationCurve;
    std::auto_ptr<fftjet::Functor1<double,double> > uncertaintyCurve;

    // Alternative method to calibrate the pileup.
    // We will fetch three lookup tables from the database:
    //
    // 1. calibration curve table
    //
    // 2. uncertainty curve table
    //
    // 3. the table that will lookup the uncertainty code
    //    given the uncalibrated pileup value
    //
    // It will be assumed that all these tables will have
    // the same record and category but different names.
    //
    std::string calibTableRecord;
    std::string calibTableCategory;
    std::string uncertaintyZonesName;
    std::string calibrationCurveName;
    std::string uncertaintyCurveName;
    bool loadCalibFromDB;
};

//
// constructors and destructor
//
FFTJetPileupEstimator::FFTJetPileupEstimator(const edm::ParameterSet& ps)
    : init_param(edm::InputTag, inputLabel),
      init_param(std::string, outputLabel),
      init_param(double, cdfvalue),
      init_param(double, ptToDensityFactor),
      init_param(unsigned, filterNumber),
      init_param(std::vector<double>, uncertaintyZones),
      init_param(std::string, calibTableRecord),
      init_param(std::string, calibTableCategory),
      init_param(std::string, uncertaintyZonesName),
      init_param(std::string, calibrationCurveName),
      init_param(std::string, uncertaintyCurveName),
      init_param(bool, loadCalibFromDB)
{
    calibrationCurve = fftjet_Function_parser(
        ps.getParameter<edm::ParameterSet>("calibrationCurve"));
    checkConfig(calibrationCurve, "bad calibration curve definition");

    uncertaintyCurve = fftjet_Function_parser(
        ps.getParameter<edm::ParameterSet>("uncertaintyCurve"));
    checkConfig(uncertaintyCurve, "bad uncertainty curve definition");

    produces<reco::FFTJetPileupSummary>(outputLabel);
}


FFTJetPileupEstimator::~FFTJetPileupEstimator()
{
}

//
// member functions
//

// ------------ method called to for each event  ------------
void FFTJetPileupEstimator::produce(edm::Event& iEvent,
                                    const edm::EventSetup& iSetup)
{
    edm::Handle<reco::DiscretizedEnergyFlow> input;
    iEvent.getByLabel(inputLabel, input);

    const reco::DiscretizedEnergyFlow& h(*input);
    const unsigned nScales = h.nEtaBins();
    const unsigned nCdfvalues = h.nPhiBins();

    const unsigned fixedCdfvalueBin = static_cast<unsigned>(
        std::floor(cdfvalue*nCdfvalues));
    if (fixedCdfvalueBin >= nCdfvalues)
    {
        throw cms::Exception("FFTJetBadConfig") 
            << "Bad cdf value" << std::endl;
    }
    if (filterNumber >= nScales)
    {
        throw cms::Exception("FFTJetBadConfig") 
            << "Bad filter number" << std::endl;
    }

    // Simple fixed-point pile-up estimate
    const double curve = h.data()[filterNumber*nCdfvalues + fixedCdfvalueBin];

    std::auto_ptr<reco::FFTJetPileupSummary> summary;
    if (loadCalibFromDB)
        summary = calibrateFromDB(curve, iSetup);
    else
        summary = calibrateFromConfig(curve);
    iEvent.put(summary, outputLabel);
}


void FFTJetPileupEstimator::beginJob()
{
}


void FFTJetPileupEstimator::endJob()
{
}


std::auto_ptr<reco::FFTJetPileupSummary>
FFTJetPileupEstimator::calibrateFromConfig(const double curve) const
{
    const double pileupRho = ptToDensityFactor*(*calibrationCurve)(curve);
    const double rhoUncert = ptToDensityFactor*(*uncertaintyCurve)(curve);

    // Determine the uncertainty zone of the estimate. The "curve"
    // has to be above or equal to uncertaintyZones[i]  but below
    // uncertaintyZones[i + 1] (the second condition is also satisfied
    // by i == uncertaintyZones.size() - 1). Of course, it is assumed
    // that the vector of zones is configured appropriately -- the zone
    // boundaries must be presented in the increasing order.
    int uncertaintyCode = -1;
    if (!uncertaintyZones.empty())
    {
        const unsigned nZones = uncertaintyZones.size();
        for (unsigned i = 0; i < nZones; ++i)
            if (curve >= uncertaintyZones[i])
            {
                if (i == nZones - 1U)
                {
                    uncertaintyCode = i;
                    break;
                }
                else if (curve < uncertaintyZones[i + 1])
                {
                    uncertaintyCode = i;
                    break;
                }
            }
    }

    return std::auto_ptr<reco::FFTJetPileupSummary>(
        new reco::FFTJetPileupSummary(curve, pileupRho,
                                      rhoUncert, uncertaintyCode));
}


std::auto_ptr<reco::FFTJetPileupSummary>
FFTJetPileupEstimator::calibrateFromDB(
    const double curve, const edm::EventSetup& iSetup) const
{
    edm::ESHandle<FFTJetLookupTableSequence> h;
    StaticFFTJetLookupTableSequenceLoader::instance().load(
        iSetup, calibTableRecord, h);
    boost::shared_ptr<npstat::StorableMultivariateFunctor> uz =
        (*h)[calibTableCategory][uncertaintyZonesName];
    boost::shared_ptr<npstat::StorableMultivariateFunctor> cc =
        (*h)[calibTableCategory][calibrationCurveName];
    boost::shared_ptr<npstat::StorableMultivariateFunctor> uc =
        (*h)[calibTableCategory][uncertaintyCurveName];

    const double pileupRho = ptToDensityFactor*(*cc)(&curve, 1U);
    const double rhoUncert = ptToDensityFactor*(*uc)(&curve, 1U);
    const int uncertaintyCode = round((*uz)(&curve, 1U));

    return std::auto_ptr<reco::FFTJetPileupSummary>(
        new reco::FFTJetPileupSummary(curve, pileupRho,
                                      rhoUncert, uncertaintyCode));
}


//define this as a plug-in
DEFINE_FWK_MODULE(FFTJetPileupEstimator);