EtaDependentPileup.cc

Go to the documentation of this file.

#include <fstream>

#include "RecoJets/FFTJetAlgorithms/interface/EtaDependentPileup.h"

namespace fftjetcms {
  EtaDependentPileup::EtaDependentPileup(const fftjet::LinearInterpolator2d& interp,
                                         const double inputRhoFactor,
                                         const double outputRhoFactor)
      : interp_(interp),
        inputRhoFactor_(inputRhoFactor),
        outputRhoFactor_(outputRhoFactor),
        etaMin_(interp.xMin()),
        etaMax_(interp.xMax()),
        rhoMin_(interp.yMin()),
        rhoMax_(interp.yMax()),
        rhoStep_((rhoMax_ - rhoMin_) / interp.ny()) {
    const double etaStep = (etaMax_ - etaMin_) / interp.nx();
    etaMin_ += etaStep * 0.5;
    etaMax_ -= etaStep * 0.5;
    rhoMin_ += rhoStep_ * 0.5;
    rhoMax_ -= rhoStep_ * 0.5;
    assert(etaMin_ < etaMax_);
    assert(rhoMin_ < rhoMax_);
  }

  double EtaDependentPileup::operator()(double eta, double /* phi */, const reco::FFTJetPileupSummary& summary) const {
    double value = 0.0;
    const double rho = summary.pileupRho() * inputRhoFactor_;
    if (eta < etaMin_)
      eta = etaMin_;
    if (eta > etaMax_)
      eta = etaMax_;
    if (rho >= rhoMin_ && rho <= rhoMax_)
      value = interp_(eta, rho);
    else {
      double x0, x1;
      if (rho < rhoMin_) {
        x0 = rhoMin_;
        x1 = rhoMin_ + rhoStep_ * 0.5;
      } else {
        x0 = rhoMax_;
        x1 = rhoMax_ - rhoStep_ * 0.5;
      }
      const double z0 = interp_(eta, x0);
      const double z1 = interp_(eta, x1);
      value = z0 + (z1 - z0) * ((rho - x0) / (x1 - x0));
    }
    return (value > 0.0 ? value : 0.0) * outputRhoFactor_;
  }
}  // namespace fftjetcms