CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_10/src/RecoJets/FFTJetProducers/src/FFTJetParameterParser.cc

Go to the documentation of this file.
00001 #include <cstdlib>
00002 #include <cstring>
00003 #include <string>
00004 #include <fstream>
00005 
00006 #include "RecoJets/FFTJetProducers/interface/FFTJetParameterParser.h"
00007 
00008 #include "RecoJets/FFTJetAlgorithms/interface/JetConvergenceDistance.h"
00009 #include "RecoJets/FFTJetAlgorithms/interface/ScaleCalculators.h"
00010 #include "RecoJets/FFTJetAlgorithms/interface/EtaAndPtDependentPeakSelector.h"
00011 #include "RecoJets/FFTJetAlgorithms/interface/EtaDependentPileup.h"
00012 #include "RecoJets/FFTJetAlgorithms/interface/PileupGrid2d.h"
00013 #include "RecoJets/FFTJetAlgorithms/interface/JetAbsEta.h"
00014 
00015 #include "FWCore/Utilities/interface/Exception.h"
00016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00017 
00018 #include "fftjet/PeakSelectors.hh"
00019 #include "fftjet/Kernel2dFactory.hh"
00020 #include "fftjet/GaussianNoiseMembershipFcn.hh"
00021 #include "fftjet/EquidistantSequence.hh"
00022 #include "fftjet/PeakEtaPhiDistance.hh"
00023 #include "fftjet/PeakEtaDependentDistance.hh"
00024 #include "fftjet/JetProperty.hh"
00025 #include "fftjet/InterpolatedMembershipFcn.hh"
00026 #include "fftjet/CompositeKernel.hh"
00027 #include "fftjet/InterpolatedKernel.hh"
00028 #include "fftjet/InterpolatedKernel3d.hh"
00029 #include "fftjet/MagneticSmearingKernel.hh"
00030 
00031 #define make_param(type, varname) const \
00032     type & varname (ps.getParameter< type >( #varname ))
00033 
00034 using namespace fftjetcms;
00035 
00036 typedef fftjet::RecombinedJet<VectorLike> RecoFFTJet;
00037 
00038 static bool parse_peak_member_function(const char* fname,
00039     fftjet::JetProperty<fftjet::Peak>::JetMemberFunction *f)
00040 {
00041     if (strcmp(fname, "eta") == 0)
00042         *f = &fftjet::Peak::eta;
00043     else if (strcmp(fname, "phi") == 0)
00044         *f = &fftjet::Peak::phi;
00045     else if (strcmp(fname, "magnitude") == 0)
00046         *f = &fftjet::Peak::magnitude;
00047     else if (strcmp(fname, "driftSpeed") == 0)
00048         *f = &fftjet::Peak::driftSpeed;
00049     else if (strcmp(fname, "magSpeed") == 0)
00050         *f = &fftjet::Peak::magSpeed;
00051     else if (strcmp(fname, "lifetime") == 0)
00052         *f = &fftjet::Peak::lifetime;
00053     else if (strcmp(fname, "scale") == 0)
00054         *f = &fftjet::Peak::scale;
00055     else if (strcmp(fname, "nearestNeighborDistance") == 0)
00056         *f = &fftjet::Peak::nearestNeighborDistance;
00057     else if (strcmp(fname, "membershipFactor") == 0)
00058         *f = &fftjet::Peak::membershipFactor;
00059     else if (strcmp(fname, "recoScale") == 0)
00060         *f = &fftjet::Peak::recoScale;
00061     else if (strcmp(fname, "recoScaleRatio") == 0)
00062         *f = &fftjet::Peak::recoScaleRatio;
00063     else if (strcmp(fname, "laplacian") == 0)
00064         *f = &fftjet::Peak::laplacian;
00065     else if (strcmp(fname, "hessianDeterminant") == 0)
00066         *f = &fftjet::Peak::hessianDeterminant;
00067     else if (strcmp(fname, "clusterRadius") == 0)
00068         *f = &fftjet::Peak::clusterRadius;
00069     else if (strcmp(fname, "clusterSeparation") == 0)
00070         *f = &fftjet::Peak::clusterSeparation;
00071     else
00072     {
00073         return false;
00074     }
00075     return true;
00076 }
00077 
00078 static bool parse_jet_member_function(const char* fname,
00079     fftjet::JetProperty<RecoFFTJet>::JetMemberFunction *f)
00080 {
00081     if (strcmp(fname, "peakEta") == 0)
00082         *f = &RecoFFTJet::peakEta;
00083     else if (strcmp(fname, "peakPhi") == 0)
00084         *f = &RecoFFTJet::peakPhi;
00085     else if (strcmp(fname, "peakMagnitude") == 0)
00086         *f = &RecoFFTJet::peakMagnitude;
00087     else if (strcmp(fname, "magnitude") == 0)
00088         *f = &RecoFFTJet::magnitude;
00089     else if (strcmp(fname, "driftSpeed") == 0)
00090         *f = &RecoFFTJet::driftSpeed;
00091     else if (strcmp(fname, "magSpeed") == 0)
00092         *f = &RecoFFTJet::magSpeed;
00093     else if (strcmp(fname, "lifetime") == 0)
00094         *f = &RecoFFTJet::lifetime;
00095     else if (strcmp(fname, "scale") == 0)
00096         *f = &RecoFFTJet::scale;
00097     else if (strcmp(fname, "nearestNeighborDistance") == 0)
00098         *f = &RecoFFTJet::nearestNeighborDistance;
00099     else if (strcmp(fname, "membershipFactor") == 0)
00100         *f = &RecoFFTJet::membershipFactor;
00101     else if (strcmp(fname, "recoScale") == 0)
00102         *f = &RecoFFTJet::recoScale;
00103     else if (strcmp(fname, "recoScaleRatio") == 0)
00104         *f = &RecoFFTJet::recoScaleRatio;
00105     else if (strcmp(fname, "laplacian") == 0)
00106         *f = &RecoFFTJet::laplacian;
00107     else if (strcmp(fname, "hessianDeterminant") == 0)
00108         *f = &RecoFFTJet::hessianDeterminant;
00109     else if (strcmp(fname, "clusterRadius") == 0)
00110         *f = &RecoFFTJet::clusterRadius;
00111     else if (strcmp(fname, "clusterSeparation") == 0)
00112         *f = &RecoFFTJet::clusterSeparation;
00113     else
00114     {
00115         return false;
00116     }
00117     return true;
00118 }
00119 
00120 
00121 namespace fftjetcms {
00122 
00123 std::auto_ptr<fftjet::Grid2d<Real> >
00124 fftjet_Grid2d_parser(const edm::ParameterSet& ps)
00125 {
00126     typedef std::auto_ptr<fftjet::Grid2d<Real> > return_type;
00127     fftjet::Grid2d<Real> *g = 0;
00128 
00129     // Check if the grid should be read from file
00130     if (ps.exists("file"))
00131     {
00132         const std::string file = ps.getParameter<std::string>("file");
00133         std::ifstream in(file.c_str(),
00134                          std::ios_base::in | std::ios_base::binary);
00135         if (!in.is_open())
00136             throw cms::Exception("FFTJetBadConfig")
00137                 << "Failed to open file " << file << std::endl;
00138         g = fftjet::Grid2d<Real>::read(in);
00139         if (g == 0)
00140             throw cms::Exception("FFTJetBadConfig")
00141                 << "Failed to read file " << file << std::endl;
00142     }
00143     else
00144     {
00145         const unsigned nEtaBins = ps.getParameter<unsigned>("nEtaBins");
00146         const Real etaMin = ps.getParameter<Real>("etaMin");
00147         const Real etaMax = ps.getParameter<Real>("etaMax");
00148         const unsigned nPhiBins = ps.getParameter<unsigned>("nPhiBins");
00149         const Real phiBin0Edge = ps.getParameter<Real>("phiBin0Edge");
00150         const std::string& title = ps.getUntrackedParameter<std::string>(
00151             "title", "");
00152 
00153         if (nEtaBins == 0 || nPhiBins == 0 || etaMin >= etaMax)
00154             return return_type(NULL);
00155         
00156         g = new fftjet::Grid2d<Real>(
00157             nEtaBins,
00158             etaMin,
00159             etaMax,
00160             nPhiBins,
00161             phiBin0Edge,
00162             title.c_str()
00163         );
00164 
00165         // Check if the grid data is provided
00166         if (ps.exists("data"))
00167         {
00168             const std::vector<Real>& data = 
00169                 ps.getParameter<std::vector<Real> >("data");
00170             if (data.size() == nEtaBins*nPhiBins)
00171                 g->blockSet(&data[0], nEtaBins, nPhiBins);
00172             else
00173             {
00174                 delete g;
00175                 g = 0;
00176             }
00177         }
00178     }
00179 
00180     return return_type(g);
00181 }
00182 
00183 
00184 std::auto_ptr<fftjet::Functor1<bool,fftjet::Peak> >
00185 fftjet_PeakSelector_parser(const edm::ParameterSet& ps)
00186 {
00187     typedef std::auto_ptr<fftjet::Functor1<bool,fftjet::Peak> > return_type;
00188 
00189     const std::string peakselector_type = ps.getParameter<std::string>(
00190         "Class");
00191 
00192     if (!peakselector_type.compare("AllPeaksPass"))
00193     {
00194         return return_type(new fftjet::AllPeaksPass());
00195     }
00196 
00197     if (!peakselector_type.compare("EtaAndPtDependentPeakSelector"))
00198     {
00199         const std::string file = ps.getParameter<std::string>("file");
00200         std::ifstream in(file.c_str(),
00201                          std::ios_base::in | std::ios_base::binary);
00202         if (!in.is_open())
00203             throw cms::Exception("FFTJetBadConfig")
00204                 << "Failed to open file " << file << std::endl;
00205         EtaAndPtDependentPeakSelector* ptr =
00206             new EtaAndPtDependentPeakSelector(in);
00207         if (!ptr->isValid())
00208             throw cms::Exception("FFTJetBadConfig")
00209                 << "Failed to read file " << file << std::endl;
00210         return return_type(ptr);
00211     }
00212 
00213     if (!peakselector_type.compare("EtaAndPtLookupPeakSelector"))
00214     {
00215         make_param(unsigned, nx);
00216         make_param(unsigned, ny);
00217         make_param(double, xmin);
00218         make_param(double, xmax);
00219         make_param(double, ymin);
00220         make_param(double, ymax);
00221         make_param(std::vector<double>, data);
00222 
00223         if (xmin >= xmax || ymin >= ymax || !nx || !ny || data.size() != nx*ny)
00224             throw cms::Exception("FFTJetBadConfig")
00225                 << "Failed to configure EtaAndPtLookupPeakSelector" << std::endl;
00226 
00227         return return_type(new EtaAndPtLookupPeakSelector(
00228                                nx, xmin, xmax, ny, ymin, ymax, data));
00229     }
00230 
00231     if (!peakselector_type.compare("SimplePeakSelector"))
00232     {
00233         const double magCut = ps.getParameter<double>("magCut");
00234         const double driftSpeedCut = ps.getParameter<double>("driftSpeedCut");
00235         const double magSpeedCut = ps.getParameter<double>("magSpeedCut");
00236         const double lifeTimeCut = ps.getParameter<double>("lifeTimeCut");
00237         const double NNDCut = ps.getParameter<double>("NNDCut");
00238         const double etaCut = ps.getParameter<double>("etaCut");
00239 
00240         return return_type(new fftjet::SimplePeakSelector(
00241             magCut, driftSpeedCut ,magSpeedCut, lifeTimeCut, NNDCut, etaCut));
00242     }
00243 
00244     if (!peakselector_type.compare("ScalePowerPeakSelector"))
00245     {
00246         const double a = ps.getParameter<double>("a");
00247         const double p = ps.getParameter<double>("p");
00248         const double b = ps.getParameter<double>("b");
00249         const double etaCut = ps.getParameter<double>("etaCut");
00250 
00251         return return_type(new fftjet::ScalePowerPeakSelector(
00252                                a, p, b, etaCut));
00253     }
00254 
00255     return return_type(NULL);
00256 }
00257 
00258 
00259 std::auto_ptr<fftjet::ScaleSpaceKernel>
00260 fftjet_MembershipFunction_parser(const edm::ParameterSet& ps)
00261 {
00262     typedef std::auto_ptr<fftjet::ScaleSpaceKernel> return_type;
00263 
00264     const std::string MembershipFunction_type = ps.getParameter<std::string>(
00265         "Class");
00266 
00267     // Parse special cases first
00268     if (!MembershipFunction_type.compare("InterpolatedMembershipFcn"))
00269     {
00270         // This is a kernel defined by a 4d (sparsified) lookup table.
00271         // Here, it is simply loaded from a file using a built-in
00272         // method from fftjet. Note that the table representation
00273         // must be native binary (this will not work on platforms with
00274         // different endianity of floating point standard).
00275         const std::string file = ps.getParameter<std::string>("file");
00276         std::ifstream in(file.c_str(),
00277                          std::ios_base::in | std::ios_base::binary);
00278         if (!in.is_open())
00279             throw cms::Exception("FFTJetBadConfig")
00280                 << "Failed to open file " << file << std::endl;
00281         return return_type(fftjet::InterpolatedMembershipFcn<float>::read(in));
00282     }
00283 
00284     if (!MembershipFunction_type.compare("Composite"))
00285     {
00286         throw cms::Exception("FFTJetBadConfig")
00287             << "Parsing of CompositeKernel objects is not implemented yet"
00288             << std::endl;
00289     }
00290 
00291     if (!MembershipFunction_type.compare("MagneticSmearing"))
00292     {
00293         // This kernel represents smearing of a jet in phi
00294         // in a magnetic field. The meaning of the parameters
00295         // is explained in the comments in the MagneticSmearingKernel.hh
00296         // header file of the fftjet package.
00297         make_param(std::vector<double>, fragmentationData);
00298         make_param(double, numeratorConst);
00299         make_param(double, charge1Fraction);
00300         make_param(double, charge2Fraction);
00301         make_param(unsigned, samplesPerBin);
00302 
00303         if (fragmentationData.empty())
00304             throw cms::Exception("FFTJetBadConfig")
00305                 << "Fragmentation function data not defined for "
00306                 "MagneticSmearingKernel" << std::endl;
00307         if (samplesPerBin < 1U)
00308             throw cms::Exception("FFTJetBadConfig")
00309                 << "Bad number of samples per bin in "
00310                 "MagneticSmearingKernel" << std::endl;
00311 
00312         fftjet::LinearInterpolator1d* fragmentationFunction = 
00313             new fftjet::LinearInterpolator1d(
00314                 &fragmentationData[0], fragmentationData.size(), 0.0, 1.0);
00315 
00316         return return_type(
00317             new fftjet::MagneticSmearingKernel<fftjet::LinearInterpolator1d>(
00318                 fragmentationFunction, numeratorConst,
00319                 charge1Fraction, charge2Fraction,
00320                 samplesPerBin, true));
00321     }
00322 
00323     if (!MembershipFunction_type.compare("Interpolated"))
00324     {
00325         // This is a kernel defined by a histogram-like 2d lookup table
00326         make_param(double, sx);
00327         make_param(double, sy);
00328         make_param(int, scalePower);
00329         make_param(unsigned, nxbins);
00330         make_param(double, xmin);
00331         make_param(double, xmax);
00332         make_param(unsigned, nybins);
00333         make_param(double, ymin);
00334         make_param(double, ymax);
00335         make_param(std::vector<double>, data);
00336 
00337         if (data.size() != nxbins*nybins)
00338             throw cms::Exception("FFTJetBadConfig")
00339                 << "Bad number of data points for Interpolated kernel"
00340                 << std::endl;
00341 
00342         return return_type(new fftjet::InterpolatedKernel(
00343                                sx, sy, scalePower,
00344                                &data[0], nxbins, xmin, xmax,
00345                                nybins, ymin, ymax));
00346     }
00347 
00348     if (!MembershipFunction_type.compare("Interpolated3d"))
00349     {
00350         // This is a kernel defined by a histogram-like 3d lookup table
00351         make_param(std::vector<double>, data);
00352         make_param(std::vector<double>, scales);
00353         make_param(bool, useLogSpaceForScale);
00354         make_param(unsigned, nxbins);
00355         make_param(double, xmin);
00356         make_param(double, xmax);
00357         make_param(unsigned, nybins);
00358         make_param(double, ymin);
00359         make_param(double, ymax);
00360 
00361         if (data.size() != nxbins*nybins*scales.size())
00362             throw cms::Exception("FFTJetBadConfig")
00363                 << "Bad number of data points for Interpolated3d kernel"
00364                 << std::endl;
00365 
00366         return return_type(new fftjet::InterpolatedKernel3d(
00367                                &data[0], scales, useLogSpaceForScale,
00368                                nxbins, xmin, xmax, nybins, ymin, ymax));
00369     }
00370 
00371     // This is not a special kernel. Try one of the classes
00372     // in the kernel factory provided by FFTJet.
00373     fftjet::DefaultKernel2dFactory factory;
00374     if (factory[MembershipFunction_type] == NULL) {
00375         return return_type(NULL);
00376     }
00377 
00378     make_param(double, sx);
00379     make_param(double, sy);
00380     make_param(int, scalePower);
00381     make_param(std::vector<double>, kernelParameters);
00382 
00383     const int n_expected = factory[MembershipFunction_type]->nParameters();
00384     if (n_expected >= 0)
00385         if (static_cast<unsigned>(n_expected) != kernelParameters.size())
00386             throw cms::Exception("FFTJetBadConfig")
00387                 << "Bad number of kernel parameters" << std::endl;
00388 
00389     return std::auto_ptr<fftjet::ScaleSpaceKernel>(
00390         factory[MembershipFunction_type]->create(
00391             sx, sy, scalePower, kernelParameters));
00392 }
00393 
00394 
00395 std::auto_ptr<AbsBgFunctor>
00396 fftjet_BgFunctor_parser(const edm::ParameterSet& ps)
00397 {
00398     const std::string bg_Membership_type = ps.getParameter<std::string>(
00399         "Class");
00400 
00401     if (!bg_Membership_type.compare("GaussianNoiseMembershipFcn"))
00402     {
00403         const double minWeight = ps.getParameter<double>("minWeight");
00404         const double prior = ps.getParameter<double>("prior");
00405         return std::auto_ptr<AbsBgFunctor>(
00406             new fftjet::GaussianNoiseMembershipFcn(minWeight,prior));
00407     }
00408 
00409     return std::auto_ptr<AbsBgFunctor>(NULL);
00410 }
00411 
00412 
00413 std::auto_ptr<std::vector<double> >
00414 fftjet_ScaleSet_parser(const edm::ParameterSet& ps)
00415 {
00416     typedef std::auto_ptr<std::vector<double> > return_type;
00417 
00418     const std::string className = ps.getParameter<std::string>("Class");
00419 
00420     if (!className.compare("EquidistantInLinearSpace") ||
00421         !className.compare("EquidistantInLogSpace"))
00422     {
00423         const double minScale = ps.getParameter<double>("minScale");
00424         const double maxScale = ps.getParameter<double>("maxScale");
00425         const unsigned nScales = ps.getParameter<unsigned>("nScales");
00426 
00427         if (minScale <= 0.0 || maxScale <= 0.0 ||
00428             nScales == 0 || minScale == maxScale)
00429             return return_type(NULL);
00430 
00431         // Can't return pointers to EquidistantInLinearSpace
00432         // or EquidistantInLogSpace directly because std::vector
00433         // destructor is not virtual.
00434         if (!className.compare("EquidistantInLinearSpace"))
00435             return return_type(new std::vector<double>(
00436                                    fftjet::EquidistantInLinearSpace(
00437                                        minScale, maxScale, nScales)));
00438         else
00439             return return_type(new std::vector<double>(
00440                                    fftjet::EquidistantInLogSpace(
00441                                        minScale, maxScale, nScales)));
00442     }
00443 
00444     if (!className.compare("UserSet"))
00445     {
00446         return_type scales(new std::vector<double>(
00447             ps.getParameter<std::vector<double> >("scales")));
00448 
00449         // Verify that all scales are positive and unique
00450         const unsigned nscales = scales->size();
00451         for (unsigned i=0; i<nscales; ++i)
00452             if ((*scales)[i] <= 0.0)
00453                 return return_type(NULL);
00454 
00455         for (unsigned i=1; i<nscales; ++i)
00456             for (unsigned j=0; j<i; ++j)
00457                 if ((*scales)[i] == (*scales)[j])
00458                     return return_type(NULL);
00459 
00460         return scales;
00461     }
00462 
00463     return return_type(NULL);
00464 }
00465 
00466 
00467 std::auto_ptr<fftjet::ClusteringTreeSparsifier<fftjet::Peak,long> >
00468 fftjet_ClusteringTreeSparsifier_parser(const edm::ParameterSet& ps)
00469 {
00470     typedef std::auto_ptr<fftjet::ClusteringTreeSparsifier<fftjet::Peak,long> > return_type;
00471 
00472     const int maxLevelNumber = ps.getParameter<int>("maxLevelNumber");
00473     const unsigned filterMask = ps.getParameter<unsigned>("filterMask");
00474     const std::vector<double> userScalesV = 
00475         ps.getParameter<std::vector<double> >("userScales");
00476     const unsigned nUserScales = userScalesV.size();
00477     const double* userScales = nUserScales ? NULL : &userScalesV[0];
00478 
00479     return return_type(
00480         new fftjet::ClusteringTreeSparsifier<fftjet::Peak,long>(
00481             maxLevelNumber,
00482             filterMask,
00483             userScales,
00484             nUserScales
00485             )
00486         );
00487 }
00488 
00489 
00490 std::auto_ptr<fftjet::AbsDistanceCalculator<fftjet::Peak> >
00491 fftjet_DistanceCalculator_parser(const edm::ParameterSet& ps)
00492 {
00493     typedef std::auto_ptr<fftjet::AbsDistanceCalculator<fftjet::Peak> > return_type;
00494 
00495     const std::string calc_type = ps.getParameter<std::string>("Class");
00496 
00497     if (!calc_type.compare("PeakEtaPhiDistance"))
00498     {
00499         const double etaToPhiBandwidthRatio = ps.getParameter<double>(
00500             "etaToPhiBandwidthRatio");
00501         return return_type(new fftjet::PeakEtaPhiDistance(etaToPhiBandwidthRatio));
00502     }
00503 
00504     if (!calc_type.compare("PeakEtaDependentDistance"))
00505     {
00506         std::auto_ptr<fftjet::LinearInterpolator1d> interp = 
00507             fftjet_LinearInterpolator1d_parser(
00508                 ps.getParameter<edm::ParameterSet>("Interpolator"));
00509         const fftjet::LinearInterpolator1d* ip = interp.get();
00510         if (ip == NULL)
00511             return return_type(NULL);
00512 
00513         // Check that the interpolator is always positive
00514         const unsigned n = ip->nx();
00515         const double* data = ip->getData();
00516         for (unsigned i=0; i<n; ++i)
00517             if (data[i] <= 0.0)
00518                 return return_type(NULL);
00519         if (ip->fLow() <= 0.0 || ip->fHigh() <= 0.0)
00520             return return_type(NULL);
00521 
00522         return return_type(new fftjet::PeakEtaDependentDistance(*ip));
00523     }
00524 
00525     return return_type(NULL);
00526 }
00527 
00528 
00529 std::auto_ptr<fftjetcms::LinInterpolatedTable1D>
00530 fftjet_LinInterpolatedTable1D_parser(const edm::ParameterSet& ps)
00531 {
00532     const double xmin = ps.getParameter<double>("xmin");
00533     const double xmax = ps.getParameter<double>("xmax");
00534     const bool leftExtrapolationLinear = 
00535         ps.getParameter<bool>("leftExtrapolationLinear");
00536     const bool rightExtrapolationLinear = 
00537         ps.getParameter<bool>("rightExtrapolationLinear");
00538     const std::vector<double> data(
00539         ps.getParameter<std::vector<double> >("data"));
00540     if (data.empty())
00541         return std::auto_ptr<fftjetcms::LinInterpolatedTable1D>(NULL);
00542     else
00543         return std::auto_ptr<fftjetcms::LinInterpolatedTable1D>(
00544             new fftjetcms::LinInterpolatedTable1D(
00545                 &data[0], data.size(), xmin, xmax,
00546                 leftExtrapolationLinear, rightExtrapolationLinear));
00547 }
00548 
00549 
00550 std::auto_ptr<fftjet::LinearInterpolator1d>
00551 fftjet_LinearInterpolator1d_parser(const edm::ParameterSet& ps)
00552 {
00553     const double xmin = ps.getParameter<double>("xmin");
00554     const double xmax = ps.getParameter<double>("xmax");
00555     const double flow = ps.getParameter<double>("flow");
00556     const double fhigh = ps.getParameter<double>("fhigh");
00557     const std::vector<double> data(
00558         ps.getParameter<std::vector<double> >("data"));
00559     if (data.empty())
00560         return std::auto_ptr<fftjet::LinearInterpolator1d>(NULL);
00561     else
00562         return std::auto_ptr<fftjet::LinearInterpolator1d>(
00563             new fftjet::LinearInterpolator1d(
00564                 &data[0], data.size(), xmin, xmax, flow, fhigh));
00565 }
00566 
00567 
00568 std::auto_ptr<fftjet::LinearInterpolator2d>
00569 fftjet_LinearInterpolator2d_parser(const edm::ParameterSet& ps)
00570 {
00571     const std::string file = ps.getParameter<std::string>("file");
00572     std::ifstream in(file.c_str(),
00573                      std::ios_base::in | std::ios_base::binary);
00574     if (!in.is_open())
00575         throw cms::Exception("FFTJetBadConfig")
00576             << "Failed to open file " << file << std::endl;
00577     fftjet::LinearInterpolator2d* ip = fftjet::LinearInterpolator2d::read(in);
00578     if (!ip)
00579         throw cms::Exception("FFTJetBadConfig")
00580             << "Failed to read file " << file << std::endl;
00581     return std::auto_ptr<fftjet::LinearInterpolator2d>(ip);
00582 }
00583 
00584 
00585 std::auto_ptr<fftjet::Functor1<double,fftjet::Peak> >
00586 fftjet_PeakFunctor_parser(const edm::ParameterSet& ps)
00587 {
00588     typedef fftjet::Functor1<double,fftjet::Peak> ptr_type;
00589     typedef std::auto_ptr<ptr_type> return_type;
00590 
00591     const std::string property_type = ps.getParameter<std::string>("Class");
00592 
00593     if (!property_type.compare("Log"))
00594     {
00595         return_type wrapped = fftjet_PeakFunctor_parser(
00596             ps.getParameter<edm::ParameterSet>("function"));
00597         ptr_type* wr = wrapped.get();
00598         if (wr)
00599         {
00600             return_type result = return_type(
00601                 new fftjet::LogProperty<fftjet::Peak>(wr, true));
00602             wrapped.release();
00603             return result;
00604         }
00605     }
00606 
00607     if (!property_type.compare("PeakProperty"))
00608     {
00609         const std::string member = ps.getParameter<std::string>("member");
00610         fftjet::JetProperty<fftjet::Peak>::JetMemberFunction fcn;
00611         if (parse_peak_member_function(member.c_str(), &fcn))
00612             return return_type(
00613                 new fftjet::JetProperty<fftjet::Peak>(fcn));
00614         else
00615             return return_type(NULL);
00616     }
00617 
00618     if (!property_type.compare("MinusScaledLaplacian"))
00619     {
00620         const double sx = ps.getParameter<double>("sx");
00621         const double sy = ps.getParameter<double>("sx");
00622         return return_type(
00623             new fftjet::MinusScaledLaplacian<fftjet::Peak>(sx, sy));
00624     }
00625 
00626     if (!property_type.compare("ScaledHessianDet"))
00627     {
00628         return return_type(
00629             new fftjet::ScaledHessianDet<fftjet::Peak>());
00630     }
00631 
00632     if (!property_type.compare("ScaledMagnitude"))
00633     {
00634         return return_type(
00635             new fftjet::ScaledMagnitude<fftjet::Peak>());
00636     }
00637 
00638     if (!property_type.compare("ScaledMagnitude2"))
00639     {
00640         return return_type(
00641             new fftjet::ScaledMagnitude2<fftjet::Peak>());
00642     }
00643 
00644     if (!property_type.compare("ConstDouble"))
00645     {
00646         const double value = ps.getParameter<double>("value");
00647         return return_type(new ConstDouble<fftjet::Peak>(value));
00648     }
00649 
00650     if (!property_type.compare("ProportionalToScale"))
00651     {
00652         const double value = ps.getParameter<double>("value");
00653         return return_type(
00654             new ProportionalToScale<fftjet::Peak>(value));
00655     }
00656 
00657     if (!property_type.compare("MultiplyByConst"))
00658     {
00659         const double factor = ps.getParameter<double>("factor");
00660         return_type function = fftjet_PeakFunctor_parser(
00661             ps.getParameter<edm::ParameterSet>("function"));
00662         ptr_type* ptr = function.get();
00663         if (ptr)
00664         {
00665             return_type result = return_type(
00666                 new MultiplyByConst<fftjet::Peak>(factor, ptr, true));
00667             function.release();
00668             return result;
00669         }
00670     }
00671 
00672     if (!property_type.compare("CompositeFunctor"))
00673     {
00674         std::auto_ptr<fftjet::Functor1<double,double> > fcn1 = 
00675             fftjet_Function_parser(
00676                 ps.getParameter<edm::ParameterSet>("function1"));
00677         return_type fcn2 = fftjet_PeakFunctor_parser(
00678             ps.getParameter<edm::ParameterSet>("function2"));
00679         fftjet::Functor1<double,double>* f1 = fcn1.get();
00680         ptr_type* f2 = fcn2.get();
00681         if (f1 && f2)
00682         {
00683             return_type result = return_type(
00684                 new CompositeFunctor<fftjet::Peak>(f1, f2, true));
00685             fcn1.release();
00686             fcn2.release();
00687             return result;
00688         }
00689     }
00690 
00691     if (!property_type.compare("ProductFunctor"))
00692     {
00693         return_type fcn1 = fftjet_PeakFunctor_parser(
00694             ps.getParameter<edm::ParameterSet>("function1"));
00695         return_type fcn2 = fftjet_PeakFunctor_parser(
00696             ps.getParameter<edm::ParameterSet>("function2"));
00697         ptr_type* f1 = fcn1.get();
00698         ptr_type* f2 = fcn2.get();
00699         if (f1 && f2)
00700         {
00701             return_type result = return_type(
00702                 new ProductFunctor<fftjet::Peak>(f1, f2, true));
00703             fcn1.release();
00704             fcn2.release();
00705             return result;
00706         }
00707     }
00708 
00709     if (!property_type.compare("MagnitudeDependent"))
00710     {
00711         std::auto_ptr<fftjet::Functor1<double,double> > fcn1 = 
00712             fftjet_Function_parser(
00713                 ps.getParameter<edm::ParameterSet>("function"));
00714         fftjet::Functor1<double,double>* f1 = fcn1.get();
00715         if (f1)
00716         {
00717             return_type result = return_type(
00718                 new MagnitudeDependent<fftjet::Peak>(f1, true));
00719             fcn1.release();
00720             return result;
00721         }
00722     }
00723 
00724     if (!property_type.compare("PeakEtaDependent"))
00725     {
00726         std::auto_ptr<fftjet::Functor1<double,double> > fcn1 = 
00727             fftjet_Function_parser(
00728                 ps.getParameter<edm::ParameterSet>("function"));
00729         fftjet::Functor1<double,double>* f1 = fcn1.get();
00730         if (f1)
00731         {
00732             return_type result = return_type(new PeakEtaDependent(f1, true));
00733             fcn1.release();
00734             return result;
00735         }
00736     }
00737 
00738     return return_type(NULL);
00739 }
00740 
00741 
00742 std::auto_ptr<fftjet::Functor1<double,fftjet::RecombinedJet<VectorLike> > >
00743 fftjet_JetFunctor_parser(const edm::ParameterSet& ps)
00744 {
00745     typedef fftjet::Functor1<double,RecoFFTJet> ptr_type;
00746     typedef std::auto_ptr<ptr_type> return_type;
00747 
00748     const std::string property_type = ps.getParameter<std::string>("Class");
00749 
00750     if (!property_type.compare("Log"))
00751     {
00752         return_type wrapped = fftjet_JetFunctor_parser(
00753             ps.getParameter<edm::ParameterSet>("function"));
00754         fftjet::Functor1<double,RecoFFTJet>* wr = wrapped.get();
00755         if (wr)
00756         {
00757             return_type result = return_type(
00758                 new fftjet::LogProperty<RecoFFTJet>(wr, true));
00759             wrapped.release();
00760             return result;
00761         }
00762     }
00763 
00764     if (!property_type.compare("JetEtaDependent"))
00765     {
00766         std::auto_ptr<fftjet::Functor1<double,double> > fcn1 = 
00767             fftjet_Function_parser(
00768                 ps.getParameter<edm::ParameterSet>("function"));
00769         fftjet::Functor1<double,double>* f1 = fcn1.get();
00770         if (f1)
00771         {
00772             return_type result = return_type(new JetEtaDependent(f1, true));
00773             fcn1.release();
00774             return result;
00775         }
00776     }
00777 
00778     if (!property_type.compare("JetProperty"))
00779     {
00780         const std::string member = ps.getParameter<std::string>("member");
00781         fftjet::JetProperty<RecoFFTJet>::JetMemberFunction fcn;
00782         if (parse_jet_member_function(member.c_str(), &fcn))
00783             return return_type(
00784                 new fftjet::JetProperty<RecoFFTJet>(fcn));
00785         else
00786             return return_type(NULL);
00787     }
00788 
00789     if (!property_type.compare("ConstDouble"))
00790     {
00791         const double value = ps.getParameter<double>("value");
00792         return return_type(new ConstDouble<RecoFFTJet>(value));
00793     }
00794 
00795     if (!property_type.compare("ProportionalToScale"))
00796     {
00797         const double value = ps.getParameter<double>("value");
00798         return return_type(
00799             new ProportionalToScale<RecoFFTJet>(value));
00800     }
00801 
00802     if (!property_type.compare("MultiplyByConst"))
00803     {
00804         const double factor = ps.getParameter<double>("factor");
00805         return_type function = fftjet_JetFunctor_parser(
00806             ps.getParameter<edm::ParameterSet>("function"));
00807         ptr_type* ptr = function.get();
00808         if (ptr)
00809         {
00810             return_type result = return_type(
00811                 new MultiplyByConst<RecoFFTJet>(factor, ptr, true));
00812             function.release();
00813             return result;
00814         }
00815     }
00816 
00817     if (!property_type.compare("CompositeFunctor"))
00818     {
00819         std::auto_ptr<fftjet::Functor1<double,double> > fcn1 = 
00820             fftjet_Function_parser(
00821                 ps.getParameter<edm::ParameterSet>("function1"));
00822         return_type fcn2 = fftjet_JetFunctor_parser(
00823             ps.getParameter<edm::ParameterSet>("function2"));
00824         fftjet::Functor1<double,double>* f1 = fcn1.get();
00825         ptr_type* f2 = fcn2.get();
00826         if (f1 && f2)
00827         {
00828             return_type result = return_type(
00829                 new CompositeFunctor<RecoFFTJet>(f1, f2, true));
00830             fcn1.release();
00831             fcn2.release();
00832             return result;
00833         }
00834     }
00835 
00836     if (!property_type.compare("ProductFunctor"))
00837     {
00838         return_type fcn1 = fftjet_JetFunctor_parser(
00839             ps.getParameter<edm::ParameterSet>("function1"));
00840         return_type fcn2 = fftjet_JetFunctor_parser(
00841             ps.getParameter<edm::ParameterSet>("function2"));
00842         ptr_type* f1 = fcn1.get();
00843         ptr_type* f2 = fcn2.get();
00844         if (f1 && f2)
00845         {
00846             return_type result = return_type(
00847                 new ProductFunctor<RecoFFTJet>(f1, f2, true));
00848             fcn1.release();
00849             fcn2.release();
00850             return result;
00851         }
00852     }
00853 
00854     if (!property_type.compare("MagnitudeDependent"))
00855     {
00856         std::auto_ptr<fftjet::Functor1<double,double> > fcn1 = 
00857             fftjet_Function_parser(
00858                 ps.getParameter<edm::ParameterSet>("function"));
00859         fftjet::Functor1<double,double>* f1 = fcn1.get();
00860         if (f1)
00861         {
00862             return_type result = return_type(
00863                 new MagnitudeDependent<RecoFFTJet>(f1, true));
00864             fcn1.release();
00865             return result;
00866         }
00867     }
00868 
00869     return return_type(NULL);
00870 }
00871 
00872 
00873 std::auto_ptr<fftjet::Functor2<double,
00874                                fftjet::RecombinedJet<VectorLike>,
00875                                fftjet::RecombinedJet<VectorLike> > >
00876 fftjet_JetDistance_parser(const edm::ParameterSet& ps)
00877 {
00878     typedef std::auto_ptr<fftjet::Functor2<
00879         double,
00880         fftjet::RecombinedJet<VectorLike>,
00881         fftjet::RecombinedJet<VectorLike> > > return_type;
00882 
00883     const std::string distance_type = ps.getParameter<std::string>(
00884         "Class");
00885 
00886     if (!distance_type.compare("JetConvergenceDistance"))
00887     {
00888         make_param(double, etaToPhiBandwidthRatio);
00889         make_param(double, relativePtBandwidth);
00890 
00891         if (etaToPhiBandwidthRatio > 0.0 && relativePtBandwidth > 0.0)
00892             return return_type(new JetConvergenceDistance(
00893                 etaToPhiBandwidthRatio, relativePtBandwidth));
00894     }
00895 
00896     return return_type(NULL);
00897 }
00898 
00899 
00900 std::auto_ptr<fftjet::Functor1<double,double> >
00901 fftjet_Function_parser(const edm::ParameterSet& ps)
00902 {
00903     typedef std::auto_ptr<fftjet::Functor1<double,double> > return_type;
00904 
00905     const std::string fcn_type = ps.getParameter<std::string>("Class");
00906 
00907     if (!fcn_type.compare("LinearInterpolator1d"))
00908     {
00909         std::auto_ptr<fftjet::LinearInterpolator1d> p = 
00910             fftjet_LinearInterpolator1d_parser(ps);
00911         fftjet::LinearInterpolator1d* ptr = p.get();
00912         if (ptr)
00913         {
00914             p.release();
00915             return return_type(ptr);
00916         }
00917     }
00918 
00919     if (!fcn_type.compare("LinInterpolatedTable1D"))
00920     {
00921         std::auto_ptr<fftjetcms::LinInterpolatedTable1D> p = 
00922             fftjet_LinInterpolatedTable1D_parser(ps);
00923         fftjetcms::LinInterpolatedTable1D* ptr = p.get();
00924         if (ptr)
00925         {
00926             p.release();
00927             return return_type(ptr);
00928         }
00929     }
00930 
00931     if (!fcn_type.compare("Polynomial"))
00932     {
00933         std::vector<double> coeffs;
00934         for (unsigned i=0; ; ++i)
00935         {
00936             std::ostringstream s;
00937             s << 'c' << i;
00938             if (ps.exists(s.str()))
00939                 coeffs.push_back(ps.getParameter<double>(s.str()));
00940             else
00941                 break;
00942         }
00943         return return_type(new Polynomial(coeffs));
00944     }
00945 
00946     return return_type(NULL);
00947 }
00948 
00949 
00950 std::auto_ptr<AbsPileupCalculator>
00951 fftjet_PileupCalculator_parser(const edm::ParameterSet& ps)
00952 {
00953     typedef std::auto_ptr<AbsPileupCalculator> return_type;
00954 
00955     const std::string fcn_type = ps.getParameter<std::string>("Class");
00956 
00957     if (!fcn_type.compare("EtaDependentPileup"))
00958     {
00959         std::auto_ptr<fftjet::LinearInterpolator2d> interp = 
00960             fftjet_LinearInterpolator2d_parser(
00961                 ps.getParameter<edm::ParameterSet>("Interpolator2d"));
00962         const double inputRhoFactor = ps.getParameter<double>("inputRhoFactor");
00963         const double outputRhoFactor = ps.getParameter<double>("outputRhoFactor");
00964 
00965         const fftjet::LinearInterpolator2d* ip = interp.get();
00966         if (ip)
00967             return return_type(new EtaDependentPileup(
00968                                    *ip, inputRhoFactor, outputRhoFactor));
00969         else
00970             return return_type(NULL);
00971     }
00972 
00973     if (!fcn_type.compare("PileupGrid2d"))
00974     {
00975         std::auto_ptr<fftjet::Grid2d<Real> > grid = 
00976             fftjet_Grid2d_parser(
00977                 ps.getParameter<edm::ParameterSet>("Grid2d"));
00978         const double rhoFactor = ps.getParameter<double>("rhoFactor");
00979 
00980         const fftjet::Grid2d<Real>* g = grid.get();
00981         if (g)
00982             return return_type(new PileupGrid2d(*g, rhoFactor));
00983         else
00984             return return_type(NULL);
00985     }
00986 
00987     return return_type(NULL);
00988 }
00989 
00990 
00991 std::auto_ptr<fftjet::JetMagnitudeMapper2d <fftjet::Peak> >
00992 fftjet_PeakMagnitudeMapper2d_parser (const edm::ParameterSet& ps)
00993 {
00994     std::auto_ptr<fftjet::LinearInterpolator2d> responseCurve =
00995         fftjet_LinearInterpolator2d_parser(
00996             ps.getParameter<edm::ParameterSet>("responseCurve"));
00997 
00998     const double minPredictor = ps.getParameter<double>("minPredictor");
00999     const double maxPredictor = ps.getParameter<double>("maxPredictor");
01000     const unsigned nPredPoints = ps.getParameter<unsigned>("nPredPoints");
01001     const double maxMagnitude = ps.getParameter<double>("maxMagnitude");
01002     const unsigned nMagPoints = ps.getParameter<unsigned>("nMagPoints");
01003 
01004     return (std::auto_ptr<fftjet::JetMagnitudeMapper2d<fftjet::Peak> >
01005              (new fftjet::JetMagnitudeMapper2d<fftjet::Peak>(
01006                   *responseCurve,
01007                   new fftjetcms::PeakAbsEta<fftjet::Peak>(),
01008                   true,minPredictor,maxPredictor,nPredPoints,
01009                   maxMagnitude,nMagPoints)));
01010 }
01011 
01012 
01013 std::auto_ptr<fftjet::JetMagnitudeMapper2d <fftjet::RecombinedJet<VectorLike> > >
01014 fftjet_JetMagnitudeMapper2d_parser (const edm::ParameterSet& ps)
01015 {
01016     std::auto_ptr<fftjet::LinearInterpolator2d> responseCurve =
01017         fftjet_LinearInterpolator2d_parser(
01018             ps.getParameter<edm::ParameterSet>("responseCurve"));
01019 
01020     const double minPredictor = ps.getParameter<double>("minPredictor");
01021     const double maxPredictor = ps.getParameter<double>("maxPredictor");
01022     const unsigned nPredPoints = ps.getParameter<unsigned>("nPredPoints");
01023     const double maxMagnitude = ps.getParameter<double>("maxMagnitude");
01024     const unsigned nMagPoints = ps.getParameter<unsigned>("nMagPoints");
01025 
01026     return (std::auto_ptr<fftjet::JetMagnitudeMapper2d<RecoFFTJet> >
01027             (new fftjet::JetMagnitudeMapper2d<RecoFFTJet>(
01028                  *responseCurve,
01029                  new fftjetcms::JetAbsEta<RecoFFTJet>(),
01030                  true,minPredictor,maxPredictor,nPredPoints,
01031                  maxMagnitude,nMagPoints)));
01032 }
01033 
01034 }