FFTJetProducers/src/FFTJetParameterParser.cc

0001 #include <cstdlib>
0002 #include <cstring>
0003 #include <memory>
0004
0005 #include <string>
0006 #include <fstream>
0007
0008 #include "RecoJets/FFTJetProducers/interface/FFTJetParameterParser.h"
0009
0010 #include "RecoJets/FFTJetAlgorithms/interface/JetConvergenceDistance.h"
0011 #include "RecoJets/FFTJetAlgorithms/interface/ScaleCalculators.h"
0012 #include "RecoJets/FFTJetAlgorithms/interface/EtaAndPtDependentPeakSelector.h"
0013 #include "RecoJets/FFTJetAlgorithms/interface/EtaDependentPileup.h"
0014 #include "RecoJets/FFTJetAlgorithms/interface/PileupGrid2d.h"
0015 #include "RecoJets/FFTJetAlgorithms/interface/JetAbsEta.h"
0016
0017 #include "FWCore/Utilities/interface/Exception.h"
0018 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0019
0020 #include "fftjet/PeakSelectors.hh"
0021 #include "fftjet/Kernel2dFactory.hh"
0022 #include "fftjet/GaussianNoiseMembershipFcn.hh"
0023 #include "fftjet/EquidistantSequence.hh"
0024 #include "fftjet/PeakEtaPhiDistance.hh"
0025 #include "fftjet/PeakEtaDependentDistance.hh"
0026 #include "fftjet/JetProperty.hh"
0027 #include "fftjet/InterpolatedMembershipFcn.hh"
0028 #include "fftjet/CompositeKernel.hh"
0029 #include "fftjet/InterpolatedKernel.hh"
0030 #include "fftjet/InterpolatedKernel3d.hh"
0031 #include "fftjet/MagneticSmearingKernel.hh"
0032
0033 #define make_param(type, varname) const type& varname(ps.getParameter<type>(#varname))
0034
0035 using namespace fftjetcms;
0036
0037 typedef fftjet::RecombinedJet<VectorLike> RecoFFTJet;
0038
0039 static bool parse_peak_member_function(const char* fname, fftjet::JetProperty<fftjet::Peak>::JetMemberFunction* f) {
0040   if (strcmp(fname, "eta") == 0)
0041     *f = &fftjet::Peak::eta;
0042   else if (strcmp(fname, "phi") == 0)
0043     *f = &fftjet::Peak::phi;
0044   else if (strcmp(fname, "magnitude") == 0)
0045     *f = &fftjet::Peak::magnitude;
0046   else if (strcmp(fname, "driftSpeed") == 0)
0047     *f = &fftjet::Peak::driftSpeed;
0048   else if (strcmp(fname, "magSpeed") == 0)
0049     *f = &fftjet::Peak::magSpeed;
0050   else if (strcmp(fname, "lifetime") == 0)
0051     *f = &fftjet::Peak::lifetime;
0052   else if (strcmp(fname, "mergeTime") == 0)
0053     *f = &fftjet::Peak::mergeTime;
0054   else if (strcmp(fname, "splitTime") == 0)
0055     *f = &fftjet::Peak::splitTime;
0056   else if (strcmp(fname, "scale") == 0)
0057     *f = &fftjet::Peak::scale;
0058   else if (strcmp(fname, "nearestNeighborDistance") == 0)
0059     *f = &fftjet::Peak::nearestNeighborDistance;
0060   else if (strcmp(fname, "membershipFactor") == 0)
0061     *f = &fftjet::Peak::membershipFactor;
0062   else if (strcmp(fname, "recoScale") == 0)
0063     *f = &fftjet::Peak::recoScale;
0064   else if (strcmp(fname, "recoScaleRatio") == 0)
0065     *f = &fftjet::Peak::recoScaleRatio;
0066   else if (strcmp(fname, "laplacian") == 0)
0067     *f = &fftjet::Peak::laplacian;
0068   else if (strcmp(fname, "hessianDeterminant") == 0)
0069     *f = &fftjet::Peak::hessianDeterminant;
0070   else if (strcmp(fname, "clusterRadius") == 0)
0071     *f = &fftjet::Peak::clusterRadius;
0072   else if (strcmp(fname, "clusterSeparation") == 0)
0073     *f = &fftjet::Peak::clusterSeparation;
0074   else {
0075     return false;
0076   }
0077   return true;
0078 }
0079
0080 static bool parse_jet_member_function(const char* fname, fftjet::JetProperty<RecoFFTJet>::JetMemberFunction* f) {
0081   if (strcmp(fname, "peakEta") == 0)
0082     *f = &RecoFFTJet::peakEta;
0083   else if (strcmp(fname, "peakPhi") == 0)
0084     *f = &RecoFFTJet::peakPhi;
0085   else if (strcmp(fname, "peakMagnitude") == 0)
0086     *f = &RecoFFTJet::peakMagnitude;
0087   else if (strcmp(fname, "magnitude") == 0)
0088     *f = &RecoFFTJet::magnitude;
0089   else if (strcmp(fname, "driftSpeed") == 0)
0090     *f = &RecoFFTJet::driftSpeed;
0091   else if (strcmp(fname, "magSpeed") == 0)
0092     *f = &RecoFFTJet::magSpeed;
0093   else if (strcmp(fname, "lifetime") == 0)
0094     *f = &RecoFFTJet::lifetime;
0095   else if (strcmp(fname, "mergeTime") == 0)
0096     *f = &RecoFFTJet::mergeTime;
0097   else if (strcmp(fname, "splitTime") == 0)
0098     *f = &RecoFFTJet::splitTime;
0099   else if (strcmp(fname, "scale") == 0)
0100     *f = &RecoFFTJet::scale;
0101   else if (strcmp(fname, "nearestNeighborDistance") == 0)
0102     *f = &RecoFFTJet::nearestNeighborDistance;
0103   else if (strcmp(fname, "membershipFactor") == 0)
0104     *f = &RecoFFTJet::membershipFactor;
0105   else if (strcmp(fname, "recoScale") == 0)
0106     *f = &RecoFFTJet::recoScale;
0107   else if (strcmp(fname, "recoScaleRatio") == 0)
0108     *f = &RecoFFTJet::recoScaleRatio;
0109   else if (strcmp(fname, "laplacian") == 0)
0110     *f = &RecoFFTJet::laplacian;
0111   else if (strcmp(fname, "hessianDeterminant") == 0)
0112     *f = &RecoFFTJet::hessianDeterminant;
0113   else if (strcmp(fname, "clusterRadius") == 0)
0114     *f = &RecoFFTJet::clusterRadius;
0115   else if (strcmp(fname, "clusterSeparation") == 0)
0116     *f = &RecoFFTJet::clusterSeparation;
0117   else {
0118     return false;
0119   }
0120   return true;
0121 }
0122
0123 namespace fftjetcms {
0124
0125   std::unique_ptr<fftjet::Grid2d<Real>> fftjet_Grid2d_parser(const edm::ParameterSet& ps) {
0126     typedef std::unique_ptr<fftjet::Grid2d<Real>> return_type;
0127     fftjet::Grid2d<Real>* g = nullptr;
0128
0129     // Check if the grid should be read from file
0130     if (ps.exists("file")) {
0131       const std::string file = ps.getParameter<std::string>("file");
0132       std::ifstream in(file.c_str(), std::ios_base::in | std::ios_base::binary);
0133       if (!in.is_open())
0134         throw cms::Exception("FFTJetBadConfig") << "Failed to open file " << file << std::endl;
0135       g = fftjet::Grid2d<Real>::read(in);
0136       if (g == nullptr)
0137         throw cms::Exception("FFTJetBadConfig") << "Failed to read file " << file << std::endl;
0138     } else {
0139       const unsigned nEtaBins = ps.getParameter<unsigned>("nEtaBins");
0140       const Real etaMin = ps.getParameter<Real>("etaMin");
0141       const Real etaMax = ps.getParameter<Real>("etaMax");
0142       const unsigned nPhiBins = ps.getParameter<unsigned>("nPhiBins");
0143       const Real phiBin0Edge = ps.getParameter<Real>("phiBin0Edge");
0144       const std::string& title = ps.getUntrackedParameter<std::string>("title", "");
0145
0146       if (nEtaBins == 0 || nPhiBins == 0 || etaMin >= etaMax)
0147         return return_type(nullptr);
0148
0149       g = new fftjet::Grid2d<Real>(nEtaBins, etaMin, etaMax, nPhiBins, phiBin0Edge, title.c_str());
0150
0151       // Check if the grid data is provided
0152       if (ps.exists("data")) {
0153         const std::vector<Real>& data = ps.getParameter<std::vector<Real>>("data");
0154         if (data.size() == nEtaBins * nPhiBins)
0155           g->blockSet(&data[0], nEtaBins, nPhiBins);
0156         else {
0157           delete g;
0158           g = nullptr;
0159         }
0160       }
0161     }
0162
0163     return return_type(g);
0164   }
0165
0166   std::unique_ptr<fftjet::Functor1<bool, fftjet::Peak>> fftjet_PeakSelector_parser(const edm::ParameterSet& ps) {
0167     typedef std::unique_ptr<fftjet::Functor1<bool, fftjet::Peak>> return_type;
0168
0169     const std::string peakselector_type = ps.getParameter<std::string>("Class");
0170
0171     if (!peakselector_type.compare("AllPeaksPass")) {
0172       return return_type(new fftjet::AllPeaksPass());
0173     }
0174
0175     if (!peakselector_type.compare("EtaAndPtDependentPeakSelector")) {
0176       const std::string file = ps.getParameter<std::string>("file");
0177       std::ifstream in(file.c_str(), std::ios_base::in | std::ios_base::binary);
0178       if (!in.is_open())
0179         throw cms::Exception("FFTJetBadConfig") << "Failed to open file " << file << std::endl;
0180       EtaAndPtDependentPeakSelector* ptr = new EtaAndPtDependentPeakSelector(in);
0181       if (!ptr->isValid())
0182         throw cms::Exception("FFTJetBadConfig") << "Failed to read file " << file << std::endl;
0183       return return_type(ptr);
0184     }
0185
0186     if (!peakselector_type.compare("EtaAndPtLookupPeakSelector")) {
0187       make_param(unsigned, nx);
0188       make_param(unsigned, ny);
0189       make_param(double, xmin);
0190       make_param(double, xmax);
0191       make_param(double, ymin);
0192       make_param(double, ymax);
0193       make_param(std::vector<double>, data);
0194
0195       if (xmin >= xmax || ymin >= ymax || !nx || !ny || data.size() != nx * ny)
0196         throw cms::Exception("FFTJetBadConfig") << "Failed to configure EtaAndPtLookupPeakSelector" << std::endl;
0197
0198       return return_type(new EtaAndPtLookupPeakSelector(nx, xmin, xmax, ny, ymin, ymax, data));
0199     }
0200
0201     if (!peakselector_type.compare("SimplePeakSelector")) {
0202       const double magCut = ps.getParameter<double>("magCut");
0203       const double driftSpeedCut = ps.getParameter<double>("driftSpeedCut");
0204       const double magSpeedCut = ps.getParameter<double>("magSpeedCut");
0205       const double lifeTimeCut = ps.getParameter<double>("lifeTimeCut");
0206       const double NNDCut = ps.getParameter<double>("NNDCut");
0207       const double etaCut = ps.getParameter<double>("etaCut");
0208       const double splitTimeCut = ps.getParameter<double>("splitTimeCut");
0209       const double mergeTimeCut = ps.getParameter<double>("mergeTimeCut");
0210
0211       return return_type(new fftjet::SimplePeakSelector(
0212           magCut, driftSpeedCut, magSpeedCut, lifeTimeCut, NNDCut, etaCut, splitTimeCut, mergeTimeCut));
0213     }
0214
0215     if (!peakselector_type.compare("ScalePowerPeakSelector")) {
0216       const double a = ps.getParameter<double>("a");
0217       const double p = ps.getParameter<double>("p");
0218       const double b = ps.getParameter<double>("b");
0219       const double etaCut = ps.getParameter<double>("etaCut");
0220
0221       return return_type(new fftjet::ScalePowerPeakSelector(a, p, b, etaCut));
0222     }
0223
0224     return return_type(nullptr);
0225   }
0226
0227   std::unique_ptr<fftjet::ScaleSpaceKernel> fftjet_MembershipFunction_parser(const edm::ParameterSet& ps) {
0228     typedef std::unique_ptr<fftjet::ScaleSpaceKernel> return_type;
0229
0230     const std::string MembershipFunction_type = ps.getParameter<std::string>("Class");
0231
0232     // Parse special cases first
0233     if (!MembershipFunction_type.compare("InterpolatedMembershipFcn")) {
0234       // This is a kernel defined by a 4d (sparsified) lookup table.
0235       // Here, it is simply loaded from a file using a built-in
0236       // method from fftjet. Note that the table representation
0237       // must be native binary (this will not work on platforms with
0238       // different endianity of floating point standard).
0239       const std::string file = ps.getParameter<std::string>("file");
0240       std::ifstream in(file.c_str(), std::ios_base::in | std::ios_base::binary);
0241       if (!in.is_open())
0242         throw cms::Exception("FFTJetBadConfig") << "Failed to open file " << file << std::endl;
0243       return return_type(fftjet::InterpolatedMembershipFcn<float>::read(in));
0244     }
0245
0246     if (!MembershipFunction_type.compare("Composite")) {
0247       throw cms::Exception("FFTJetBadConfig")
0248           << "Parsing of CompositeKernel objects is not implemented yet" << std::endl;
0249     }
0250
0251     if (!MembershipFunction_type.compare("MagneticSmearing")) {
0252       // This kernel represents smearing of a jet in phi
0253       // in a magnetic field. The meaning of the parameters
0254       // is explained in the comments in the MagneticSmearingKernel.hh
0255       // header file of the fftjet package.
0256       make_param(std::vector<double>, fragmentationData);
0257       make_param(double, numeratorConst);
0258       make_param(double, charge1Fraction);
0259       make_param(double, charge2Fraction);
0260       make_param(unsigned, samplesPerBin);
0261
0262       if (fragmentationData.empty())
0263         throw cms::Exception("FFTJetBadConfig") << "Fragmentation function data not defined for "
0264                                                    "MagneticSmearingKernel"
0265                                                 << std::endl;
0266       if (samplesPerBin < 1U)
0267         throw cms::Exception("FFTJetBadConfig") << "Bad number of samples per bin in "
0268                                                    "MagneticSmearingKernel"
0269                                                 << std::endl;
0270
0271       fftjet::LinearInterpolator1d* fragmentationFunction =
0272           new fftjet::LinearInterpolator1d(&fragmentationData[0], fragmentationData.size(), 0.0, 1.0);
0273
0274       return return_type(new fftjet::MagneticSmearingKernel<fftjet::LinearInterpolator1d>(
0275           fragmentationFunction, numeratorConst, charge1Fraction, charge2Fraction, samplesPerBin, true));
0276     }
0277
0278     if (!MembershipFunction_type.compare("Interpolated")) {
0279       // This is a kernel defined by a histogram-like 2d lookup table
0280       make_param(double, sx);
0281       make_param(double, sy);
0282       make_param(int, scalePower);
0283       make_param(unsigned, nxbins);
0284       make_param(double, xmin);
0285       make_param(double, xmax);
0286       make_param(unsigned, nybins);
0287       make_param(double, ymin);
0288       make_param(double, ymax);
0289       make_param(std::vector<double>, data);
0290
0291       if (data.size() != nxbins * nybins)
0292         throw cms::Exception("FFTJetBadConfig") << "Bad number of data points for Interpolated kernel" << std::endl;
0293
0294       return return_type(
0295           new fftjet::InterpolatedKernel(sx, sy, scalePower, &data[0], nxbins, xmin, xmax, nybins, ymin, ymax));
0296     }
0297
0298     if (!MembershipFunction_type.compare("Interpolated3d")) {
0299       // This is a kernel defined by a histogram-like 3d lookup table
0300       make_param(std::vector<double>, data);
0301       make_param(std::vector<double>, scales);
0302       make_param(bool, useLogSpaceForScale);
0303       make_param(unsigned, nxbins);
0304       make_param(double, xmin);
0305       make_param(double, xmax);
0306       make_param(unsigned, nybins);
0307       make_param(double, ymin);
0308       make_param(double, ymax);
0309
0310       if (data.size() != nxbins * nybins * scales.size())
0311         throw cms::Exception("FFTJetBadConfig") << "Bad number of data points for Interpolated3d kernel" << std::endl;
0312
0313       return return_type(new fftjet::InterpolatedKernel3d(
0314           &data[0], scales, useLogSpaceForScale, nxbins, xmin, xmax, nybins, ymin, ymax));
0315     }
0316
0317     // This is not a special kernel. Try one of the classes
0318     // in the kernel factory provided by FFTJet.
0319     fftjet::DefaultKernel2dFactory factory;
0320     if (factory[MembershipFunction_type] == nullptr) {
0321       return return_type(nullptr);
0322     }
0323
0324     make_param(double, sx);
0325     make_param(double, sy);
0326     make_param(int, scalePower);
0327     make_param(std::vector<double>, kernelParameters);
0328
0329     const int n_expected = factory[MembershipFunction_type]->nParameters();
0330     if (n_expected >= 0)
0331       if (static_cast<unsigned>(n_expected) != kernelParameters.size())
0332         throw cms::Exception("FFTJetBadConfig") << "Bad number of kernel parameters" << std::endl;
0333
0334     return std::unique_ptr<fftjet::ScaleSpaceKernel>(
0335         factory[MembershipFunction_type]->create(sx, sy, scalePower, kernelParameters));
0336   }
0337
0338   std::unique_ptr<AbsBgFunctor> fftjet_BgFunctor_parser(const edm::ParameterSet& ps) {
0339     const std::string bg_Membership_type = ps.getParameter<std::string>("Class");
0340
0341     if (!bg_Membership_type.compare("GaussianNoiseMembershipFcn")) {
0342       const double minWeight = ps.getParameter<double>("minWeight");
0343       const double prior = ps.getParameter<double>("prior");
0344       return std::unique_ptr<AbsBgFunctor>(new fftjet::GaussianNoiseMembershipFcn(minWeight, prior));
0345     }
0346
0347     return std::unique_ptr<AbsBgFunctor>(nullptr);
0348   }
0349
0350   std::unique_ptr<std::vector<double>> fftjet_ScaleSet_parser(const edm::ParameterSet& ps) {
0351     typedef std::unique_ptr<std::vector<double>> return_type;
0352
0353     const std::string className = ps.getParameter<std::string>("Class");
0354
0355     if (!className.compare("EquidistantInLinearSpace") || !className.compare("EquidistantInLogSpace")) {
0356       const double minScale = ps.getParameter<double>("minScale");
0357       const double maxScale = ps.getParameter<double>("maxScale");
0358       const unsigned nScales = ps.getParameter<unsigned>("nScales");
0359
0360       if (minScale <= 0.0 || maxScale <= 0.0 || nScales == 0 || minScale == maxScale)
0361         return return_type(nullptr);
0362
0363       // Can't return pointers to EquidistantInLinearSpace
0364       // or EquidistantInLogSpace directly because std::vector
0365       // destructor is not virtual.
0366       if (!className.compare("EquidistantInLinearSpace"))
0367         return std::make_unique<std::vector<double>>(fftjet::EquidistantInLinearSpace(minScale, maxScale, nScales));
0368       else
0369         return std::make_unique<std::vector<double>>(fftjet::EquidistantInLogSpace(minScale, maxScale, nScales));
0370     }
0371
0372     if (!className.compare("UserSet")) {
0373       return_type scales(new std::vector<double>(ps.getParameter<std::vector<double>>("scales")));
0374
0375       // Verify that all scales are positive and unique
0376       const unsigned nscales = scales->size();
0377       for (unsigned i = 0; i < nscales; ++i)
0378         if ((*scales)[i] <= 0.0)
0379           return return_type(nullptr);
0380
0381       for (unsigned i = 1; i < nscales; ++i)
0382         for (unsigned j = 0; j < i; ++j)
0383           if ((*scales)[i] == (*scales)[j])
0384             return return_type(nullptr);
0385
0386       return scales;
0387     }
0388
0389     return return_type(nullptr);
0390   }
0391
0392   std::unique_ptr<fftjet::ClusteringTreeSparsifier<fftjet::Peak, long>> fftjet_ClusteringTreeSparsifier_parser(
0393       const edm::ParameterSet& ps) {
0394     const int maxLevelNumber = ps.getParameter<int>("maxLevelNumber");
0395     const unsigned filterMask = ps.getParameter<unsigned>("filterMask");
0396     const std::vector<double> userScalesV = ps.getParameter<std::vector<double>>("userScales");
0397     const unsigned nUserScales = userScalesV.size();
0398     const double* userScales = nUserScales ? &userScalesV[0] : nullptr;
0399
0400     return std::make_unique<fftjet::ClusteringTreeSparsifier<fftjet::Peak, long>>(
0401         maxLevelNumber, filterMask, userScales, nUserScales);
0402   }
0403
0404   std::unique_ptr<fftjet::AbsDistanceCalculator<fftjet::Peak>> fftjet_DistanceCalculator_parser(
0405       const edm::ParameterSet& ps) {
0406     typedef std::unique_ptr<fftjet::AbsDistanceCalculator<fftjet::Peak>> return_type;
0407
0408     const std::string calc_type = ps.getParameter<std::string>("Class");
0409
0410     if (!calc_type.compare("PeakEtaPhiDistance")) {
0411       const double etaToPhiBandwidthRatio = ps.getParameter<double>("etaToPhiBandwidthRatio");
0412       return return_type(new fftjet::PeakEtaPhiDistance(etaToPhiBandwidthRatio));
0413     }
0414
0415     if (!calc_type.compare("PeakEtaDependentDistance")) {
0416       std::unique_ptr<fftjet::LinearInterpolator1d> interp =
0417           fftjet_LinearInterpolator1d_parser(ps.getParameter<edm::ParameterSet>("Interpolator"));
0418       const fftjet::LinearInterpolator1d* ip = interp.get();
0419       if (ip == nullptr)
0420         return return_type(nullptr);
0421
0422       // Check that the interpolator is always positive
0423       const unsigned n = ip->nx();
0424       const double* data = ip->getData();
0425       for (unsigned i = 0; i < n; ++i)
0426         if (data[i] <= 0.0)
0427           return return_type(nullptr);
0428       if (ip->fLow() <= 0.0 || ip->fHigh() <= 0.0)
0429         return return_type(nullptr);
0430
0431       return return_type(new fftjet::PeakEtaDependentDistance(*ip));
0432     }
0433
0434     return return_type(nullptr);
0435   }
0436
0437   std::unique_ptr<fftjetcms::LinInterpolatedTable1D> fftjet_LinInterpolatedTable1D_parser(const edm::ParameterSet& ps) {
0438     const double xmin = ps.getParameter<double>("xmin");
0439     const double xmax = ps.getParameter<double>("xmax");
0440     const bool leftExtrapolationLinear = ps.getParameter<bool>("leftExtrapolationLinear");
0441     const bool rightExtrapolationLinear = ps.getParameter<bool>("rightExtrapolationLinear");
0442     const std::vector<double> data(ps.getParameter<std::vector<double>>("data"));
0443     if (data.empty())
0444       return std::unique_ptr<fftjetcms::LinInterpolatedTable1D>(nullptr);
0445     else
0446       return std::make_unique<fftjetcms::LinInterpolatedTable1D>(
0447           &data[0], data.size(), xmin, xmax, leftExtrapolationLinear, rightExtrapolationLinear);
0448   }
0449
0450   std::unique_ptr<fftjet::LinearInterpolator1d> fftjet_LinearInterpolator1d_parser(const edm::ParameterSet& ps) {
0451     const double xmin = ps.getParameter<double>("xmin");
0452     const double xmax = ps.getParameter<double>("xmax");
0453     const double flow = ps.getParameter<double>("flow");
0454     const double fhigh = ps.getParameter<double>("fhigh");
0455     const std::vector<double> data(ps.getParameter<std::vector<double>>("data"));
0456     if (data.empty())
0457       return std::unique_ptr<fftjet::LinearInterpolator1d>(nullptr);
0458     else
0459       return std::make_unique<fftjet::LinearInterpolator1d>(&data[0], data.size(), xmin, xmax, flow, fhigh);
0460   }
0461
0462   std::unique_ptr<fftjet::LinearInterpolator2d> fftjet_LinearInterpolator2d_parser(const edm::ParameterSet& ps) {
0463     const std::string file = ps.getParameter<std::string>("file");
0464     std::ifstream in(file.c_str(), std::ios_base::in | std::ios_base::binary);
0465     if (!in.is_open())
0466       throw cms::Exception("FFTJetBadConfig") << "Failed to open file " << file << std::endl;
0467     fftjet::LinearInterpolator2d* ip = fftjet::LinearInterpolator2d::read(in);
0468     if (!ip)
0469       throw cms::Exception("FFTJetBadConfig") << "Failed to read file " << file << std::endl;
0470     return std::unique_ptr<fftjet::LinearInterpolator2d>(ip);
0471   }
0472
0473   std::unique_ptr<fftjet::Functor1<double, fftjet::Peak>> fftjet_PeakFunctor_parser(const edm::ParameterSet& ps) {
0474     typedef fftjet::Functor1<double, fftjet::Peak> ptr_type;
0475     typedef std::unique_ptr<ptr_type> return_type;
0476
0477     const std::string property_type = ps.getParameter<std::string>("Class");
0478
0479     if (!property_type.compare("Log")) {
0480       return_type wrapped = fftjet_PeakFunctor_parser(ps.getParameter<edm::ParameterSet>("function"));
0481       ptr_type* wr = wrapped.get();
0482       if (wr) {
0483         return_type result = return_type(new fftjet::LogProperty<fftjet::Peak>(wr, true));
0484         wrapped.release();
0485         return result;
0486       }
0487     }
0488
0489     if (!property_type.compare("PeakProperty")) {
0490       const std::string member = ps.getParameter<std::string>("member");
0491       fftjet::JetProperty<fftjet::Peak>::JetMemberFunction fcn;
0492       if (parse_peak_member_function(member.c_str(), &fcn))
0493         return return_type(new fftjet::JetProperty<fftjet::Peak>(fcn));
0494       else
0495         return return_type(nullptr);
0496     }
0497
0498     if (!property_type.compare("MinusScaledLaplacian")) {
0499       const double sx = ps.getParameter<double>("sx");
0500       const double sy = ps.getParameter<double>("sx");
0501       return return_type(new fftjet::MinusScaledLaplacian<fftjet::Peak>(sx, sy));
0502     }
0503
0504     if (!property_type.compare("ScaledHessianDet")) {
0505       return return_type(new fftjet::ScaledHessianDet<fftjet::Peak>());
0506     }
0507
0508     if (!property_type.compare("ScaledMagnitude")) {
0509       return return_type(new fftjet::ScaledMagnitude<fftjet::Peak>());
0510     }
0511
0512     if (!property_type.compare("ScaledMagnitude2")) {
0513       return return_type(new fftjet::ScaledMagnitude2<fftjet::Peak>());
0514     }
0515
0516     if (!property_type.compare("ConstDouble")) {
0517       const double value = ps.getParameter<double>("value");
0518       return return_type(new ConstDouble<fftjet::Peak>(value));
0519     }
0520
0521     if (!property_type.compare("ProportionalToScale")) {
0522       const double value = ps.getParameter<double>("value");
0523       return return_type(new ProportionalToScale<fftjet::Peak>(value));
0524     }
0525
0526     if (!property_type.compare("MultiplyByConst")) {
0527       const double factor = ps.getParameter<double>("factor");
0528       return_type function = fftjet_PeakFunctor_parser(ps.getParameter<edm::ParameterSet>("function"));
0529       ptr_type* ptr = function.get();
0530       if (ptr) {
0531         return_type result = return_type(new MultiplyByConst<fftjet::Peak>(factor, ptr, true));
0532         function.release();
0533         return result;
0534       }
0535     }
0536
0537     if (!property_type.compare("CompositeFunctor")) {
0538       std::unique_ptr<fftjet::Functor1<double, double>> fcn1 =
0539           fftjet_Function_parser(ps.getParameter<edm::ParameterSet>("function1"));
0540       return_type fcn2 = fftjet_PeakFunctor_parser(ps.getParameter<edm::ParameterSet>("function2"));
0541       fftjet::Functor1<double, double>* f1 = fcn1.get();
0542       ptr_type* f2 = fcn2.get();
0543       if (f1 && f2) {
0544         return_type result = return_type(new CompositeFunctor<fftjet::Peak>(f1, f2, true));
0545         fcn1.release();
0546         fcn2.release();
0547         return result;
0548       }
0549     }
0550
0551     if (!property_type.compare("ProductFunctor")) {
0552       return_type fcn1 = fftjet_PeakFunctor_parser(ps.getParameter<edm::ParameterSet>("function1"));
0553       return_type fcn2 = fftjet_PeakFunctor_parser(ps.getParameter<edm::ParameterSet>("function2"));
0554       ptr_type* f1 = fcn1.get();
0555       ptr_type* f2 = fcn2.get();
0556       if (f1 && f2) {
0557         return_type result = return_type(new ProductFunctor<fftjet::Peak>(f1, f2, true));
0558         fcn1.release();
0559         fcn2.release();
0560         return result;
0561       }
0562     }
0563
0564     if (!property_type.compare("MagnitudeDependent")) {
0565       std::unique_ptr<fftjet::Functor1<double, double>> fcn1 =
0566           fftjet_Function_parser(ps.getParameter<edm::ParameterSet>("function"));
0567       fftjet::Functor1<double, double>* f1 = fcn1.get();
0568       if (f1) {
0569         return_type result = return_type(new MagnitudeDependent<fftjet::Peak>(f1, true));
0570         fcn1.release();
0571         return result;
0572       }
0573     }
0574
0575     if (!property_type.compare("PeakEtaDependent")) {
0576       std::unique_ptr<fftjet::Functor1<double, double>> fcn1 =
0577           fftjet_Function_parser(ps.getParameter<edm::ParameterSet>("function"));
0578       fftjet::Functor1<double, double>* f1 = fcn1.get();
0579       if (f1) {
0580         return_type result = return_type(new PeakEtaDependent(f1, true));
0581         fcn1.release();
0582         return result;
0583       }
0584     }
0585
0586     return return_type(nullptr);
0587   }
0588
0589   std::unique_ptr<fftjet::Functor1<double, fftjet::RecombinedJet<VectorLike>>> fftjet_JetFunctor_parser(
0590       const edm::ParameterSet& ps) {
0591     typedef fftjet::Functor1<double, RecoFFTJet> ptr_type;
0592     typedef std::unique_ptr<ptr_type> return_type;
0593
0594     const std::string property_type = ps.getParameter<std::string>("Class");
0595
0596     if (!property_type.compare("Log")) {
0597       return_type wrapped = fftjet_JetFunctor_parser(ps.getParameter<edm::ParameterSet>("function"));
0598       fftjet::Functor1<double, RecoFFTJet>* wr = wrapped.get();
0599       if (wr) {
0600         return_type result = return_type(new fftjet::LogProperty<RecoFFTJet>(wr, true));
0601         wrapped.release();
0602         return result;
0603       }
0604     }
0605
0606     if (!property_type.compare("JetEtaDependent")) {
0607       std::unique_ptr<fftjet::Functor1<double, double>> fcn1 =
0608           fftjet_Function_parser(ps.getParameter<edm::ParameterSet>("function"));
0609       fftjet::Functor1<double, double>* f1 = fcn1.get();
0610       if (f1) {
0611         return_type result = return_type(new JetEtaDependent(f1, true));
0612         fcn1.release();
0613         return result;
0614       }
0615     }
0616
0617     if (!property_type.compare("JetProperty")) {
0618       const std::string member = ps.getParameter<std::string>("member");
0619       fftjet::JetProperty<RecoFFTJet>::JetMemberFunction fcn;
0620       if (parse_jet_member_function(member.c_str(), &fcn))
0621         return return_type(new fftjet::JetProperty<RecoFFTJet>(fcn));
0622       else
0623         return return_type(nullptr);
0624     }
0625
0626     if (!property_type.compare("ConstDouble")) {
0627       const double value = ps.getParameter<double>("value");
0628       return return_type(new ConstDouble<RecoFFTJet>(value));
0629     }
0630
0631     if (!property_type.compare("ProportionalToScale")) {
0632       const double value = ps.getParameter<double>("value");
0633       return return_type(new ProportionalToScale<RecoFFTJet>(value));
0634     }
0635
0636     if (!property_type.compare("MultiplyByConst")) {
0637       const double factor = ps.getParameter<double>("factor");
0638       return_type function = fftjet_JetFunctor_parser(ps.getParameter<edm::ParameterSet>("function"));
0639       ptr_type* ptr = function.get();
0640       if (ptr) {
0641         return_type result = return_type(new MultiplyByConst<RecoFFTJet>(factor, ptr, true));
0642         function.release();
0643         return result;
0644       }
0645     }
0646
0647     if (!property_type.compare("CompositeFunctor")) {
0648       std::unique_ptr<fftjet::Functor1<double, double>> fcn1 =
0649           fftjet_Function_parser(ps.getParameter<edm::ParameterSet>("function1"));
0650       return_type fcn2 = fftjet_JetFunctor_parser(ps.getParameter<edm::ParameterSet>("function2"));
0651       fftjet::Functor1<double, double>* f1 = fcn1.get();
0652       ptr_type* f2 = fcn2.get();
0653       if (f1 && f2) {
0654         return_type result = return_type(new CompositeFunctor<RecoFFTJet>(f1, f2, true));
0655         fcn1.release();
0656         fcn2.release();
0657         return result;
0658       }
0659     }
0660
0661     if (!property_type.compare("ProductFunctor")) {
0662       return_type fcn1 = fftjet_JetFunctor_parser(ps.getParameter<edm::ParameterSet>("function1"));
0663       return_type fcn2 = fftjet_JetFunctor_parser(ps.getParameter<edm::ParameterSet>("function2"));
0664       ptr_type* f1 = fcn1.get();
0665       ptr_type* f2 = fcn2.get();
0666       if (f1 && f2) {
0667         return_type result = return_type(new ProductFunctor<RecoFFTJet>(f1, f2, true));
0668         fcn1.release();
0669         fcn2.release();
0670         return result;
0671       }
0672     }
0673
0674     if (!property_type.compare("MagnitudeDependent")) {
0675       std::unique_ptr<fftjet::Functor1<double, double>> fcn1 =
0676           fftjet_Function_parser(ps.getParameter<edm::ParameterSet>("function"));
0677       fftjet::Functor1<double, double>* f1 = fcn1.get();
0678       if (f1) {
0679         return_type result = return_type(new MagnitudeDependent<RecoFFTJet>(f1, true));
0680         fcn1.release();
0681         return result;
0682       }
0683     }
0684
0685     return return_type(nullptr);
0686   }
0687
0688   std::unique_ptr<fftjet::Functor2<double, fftjet::RecombinedJet<VectorLike>, fftjet::RecombinedJet<VectorLike>>>
0689   fftjet_JetDistance_parser(const edm::ParameterSet& ps) {
0690     typedef std::unique_ptr<
0691         fftjet::Functor2<double, fftjet::RecombinedJet<VectorLike>, fftjet::RecombinedJet<VectorLike>>>
0692         return_type;
0693
0694     const std::string distance_type = ps.getParameter<std::string>("Class");
0695
0696     if (!distance_type.compare("JetConvergenceDistance")) {
0697       make_param(double, etaToPhiBandwidthRatio);
0698       make_param(double, relativePtBandwidth);
0699
0700       if (etaToPhiBandwidthRatio > 0.0 && relativePtBandwidth > 0.0)
0701         return return_type(new JetConvergenceDistance(etaToPhiBandwidthRatio, relativePtBandwidth));
0702     }
0703
0704     return return_type(nullptr);
0705   }
0706
0707   std::unique_ptr<fftjet::Functor1<double, double>> fftjet_Function_parser(const edm::ParameterSet& ps) {
0708     typedef std::unique_ptr<fftjet::Functor1<double, double>> return_type;
0709
0710     const std::string fcn_type = ps.getParameter<std::string>("Class");
0711
0712     if (!fcn_type.compare("LinearInterpolator1d")) {
0713       std::unique_ptr<fftjet::LinearInterpolator1d> p = fftjet_LinearInterpolator1d_parser(ps);
0714       fftjet::LinearInterpolator1d* ptr = p.get();
0715       if (ptr) {
0716         p.release();
0717         return return_type(ptr);
0718       }
0719     }
0720
0721     if (!fcn_type.compare("LinInterpolatedTable1D")) {
0722       std::unique_ptr<fftjetcms::LinInterpolatedTable1D> p = fftjet_LinInterpolatedTable1D_parser(ps);
0723       fftjetcms::LinInterpolatedTable1D* ptr = p.get();
0724       if (ptr) {
0725         p.release();
0726         return return_type(ptr);
0727       }
0728     }
0729
0730     if (!fcn_type.compare("Polynomial")) {
0731       std::vector<double> coeffs;
0732       for (unsigned i = 0;; ++i) {
0733         std::ostringstream s;
0734         s << 'c' << i;
0735         if (ps.exists(s.str()))
0736           coeffs.push_back(ps.getParameter<double>(s.str()));
0737         else
0738           break;
0739       }
0740       return return_type(new Polynomial(coeffs));
0741     }
0742
0743     return return_type(nullptr);
0744   }
0745
0746   std::unique_ptr<AbsPileupCalculator> fftjet_PileupCalculator_parser(const edm::ParameterSet& ps) {
0747     typedef std::unique_ptr<AbsPileupCalculator> return_type;
0748
0749     const std::string fcn_type = ps.getParameter<std::string>("Class");
0750
0751     if (!fcn_type.compare("EtaDependentPileup")) {
0752       std::unique_ptr<fftjet::LinearInterpolator2d> interp =
0753           fftjet_LinearInterpolator2d_parser(ps.getParameter<edm::ParameterSet>("Interpolator2d"));
0754       const double inputRhoFactor = ps.getParameter<double>("inputRhoFactor");
0755       const double outputRhoFactor = ps.getParameter<double>("outputRhoFactor");
0756
0757       const fftjet::LinearInterpolator2d* ip = interp.get();
0758       if (ip)
0759         return return_type(new EtaDependentPileup(*ip, inputRhoFactor, outputRhoFactor));
0760       else
0761         return return_type(nullptr);
0762     }
0763
0764     if (!fcn_type.compare("PileupGrid2d")) {
0765       std::unique_ptr<fftjet::Grid2d<Real>> grid = fftjet_Grid2d_parser(ps.getParameter<edm::ParameterSet>("Grid2d"));
0766       const double rhoFactor = ps.getParameter<double>("rhoFactor");
0767
0768       const fftjet::Grid2d<Real>* g = grid.get();
0769       if (g)
0770         return return_type(new PileupGrid2d(*g, rhoFactor));
0771       else
0772         return return_type(nullptr);
0773     }
0774
0775     return return_type(nullptr);
0776   }
0777
0778   std::unique_ptr<fftjet::JetMagnitudeMapper2d<fftjet::Peak>> fftjet_PeakMagnitudeMapper2d_parser(
0779       const edm::ParameterSet& ps) {
0780     std::unique_ptr<fftjet::LinearInterpolator2d> responseCurve =
0781         fftjet_LinearInterpolator2d_parser(ps.getParameter<edm::ParameterSet>("responseCurve"));
0782
0783     const double minPredictor = ps.getParameter<double>("minPredictor");
0784     const double maxPredictor = ps.getParameter<double>("maxPredictor");
0785     const unsigned nPredPoints = ps.getParameter<unsigned>("nPredPoints");
0786     const double maxMagnitude = ps.getParameter<double>("maxMagnitude");
0787     const unsigned nMagPoints = ps.getParameter<unsigned>("nMagPoints");
0788
0789     return (std::make_unique<fftjet::JetMagnitudeMapper2d<fftjet::Peak>>(*responseCurve,
0790                                                                          new fftjetcms::PeakAbsEta<fftjet::Peak>(),
0791                                                                          true,
0792                                                                          minPredictor,
0793                                                                          maxPredictor,
0794                                                                          nPredPoints,
0795                                                                          maxMagnitude,
0796                                                                          nMagPoints));
0797   }
0798
0799   std::unique_ptr<fftjet::JetMagnitudeMapper2d<fftjet::RecombinedJet<VectorLike>>> fftjet_JetMagnitudeMapper2d_parser(
0800       const edm::ParameterSet& ps) {
0801     std::unique_ptr<fftjet::LinearInterpolator2d> responseCurve =
0802         fftjet_LinearInterpolator2d_parser(ps.getParameter<edm::ParameterSet>("responseCurve"));
0803
0804     const double minPredictor = ps.getParameter<double>("minPredictor");
0805     const double maxPredictor = ps.getParameter<double>("maxPredictor");
0806     const unsigned nPredPoints = ps.getParameter<unsigned>("nPredPoints");
0807     const double maxMagnitude = ps.getParameter<double>("maxMagnitude");
0808     const unsigned nMagPoints = ps.getParameter<unsigned>("nMagPoints");
0809
0810     return (std::make_unique<fftjet::JetMagnitudeMapper2d<RecoFFTJet>>(*responseCurve,
0811                                                                        new fftjetcms::JetAbsEta<RecoFFTJet>(),
0812                                                                        true,
0813                                                                        minPredictor,
0814                                                                        maxPredictor,
0815                                                                        nPredPoints,
0816                                                                        maxMagnitude,
0817                                                                        nMagPoints));
0818   }
0819
0820 }  // namespace fftjetcms