FFTJetAlgorithms/src/EtaDependentPileup.cc

0001 #include <fstream>
0002
0003 #include "RecoJets/FFTJetAlgorithms/interface/EtaDependentPileup.h"
0004
0005 namespace fftjetcms {
0006   EtaDependentPileup::EtaDependentPileup(const fftjet::LinearInterpolator2d& interp,
0007                                          const double inputRhoFactor,
0008                                          const double outputRhoFactor)
0009       : interp_(interp),
0010         inputRhoFactor_(inputRhoFactor),
0011         outputRhoFactor_(outputRhoFactor),
0012         etaMin_(interp.xMin()),
0013         etaMax_(interp.xMax()),
0014         rhoMin_(interp.yMin()),
0015         rhoMax_(interp.yMax()),
0016         rhoStep_((rhoMax_ - rhoMin_) / interp.ny()) {
0017     const double etaStep = (etaMax_ - etaMin_) / interp.nx();
0018     etaMin_ += etaStep * 0.5;
0019     etaMax_ -= etaStep * 0.5;
0020     rhoMin_ += rhoStep_ * 0.5;
0021     rhoMax_ -= rhoStep_ * 0.5;
0022     assert(etaMin_ < etaMax_);
0023     assert(rhoMin_ < rhoMax_);
0024   }
0025
0026   double EtaDependentPileup::operator()(double eta, double /* phi */, const reco::FFTJetPileupSummary& summary) const {
0027     double value = 0.0;
0028     const double rho = summary.pileupRho() * inputRhoFactor_;
0029     if (eta < etaMin_)
0030       eta = etaMin_;
0031     if (eta > etaMax_)
0032       eta = etaMax_;
0033     if (rho >= rhoMin_ && rho <= rhoMax_)
0034       value = interp_(eta, rho);
0035     else {
0036       double x0, x1;
0037       if (rho < rhoMin_) {
0038         x0 = rhoMin_;
0039         x1 = rhoMin_ + rhoStep_ * 0.5;
0040       } else {
0041         x0 = rhoMax_;
0042         x1 = rhoMax_ - rhoStep_ * 0.5;
0043       }
0044       const double z0 = interp_(eta, x0);
0045       const double z1 = interp_(eta, x1);
0046       value = z0 + (z1 - z0) * ((rho - x0) / (x1 - x0));
0047     }
0048     return (value > 0.0 ? value : 0.0) * outputRhoFactor_;
0049   }
0050 }  // namespace fftjetcms