File indexing completed on 2024-04-06 12:25:22
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 , 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 }