Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-05-21 03:43:06 

0001 // -*- C++ -*-
0002 //
0003 // Package:    FFTJetProducers
0004 // Class:      FFTJetProducer
0005 //
0006 /**\class FFTJetProducer FFTJetProducer.cc RecoJets/FFTJetProducers/plugins/FFTJetProducer.cc
0007 
0008  Description: makes jets using FFTJet clustering tree
0009 
0010  Implementation:
0011      [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  Igor Volobouev
0015 //         Created:  Sun Jun 20 14:32:36 CDT 2010
0016 //
0017 //
0018 
0019 #include <algorithm>
0020 #include <fstream>
0021 #include <functional>
0022 #include <iostream>
0023 #include <memory>
0024 
0025 #include "fftjet/peakEtLifetime.hh"
0026 
0027 // Header for this class
0028 #include "RecoJets/FFTJetProducers/plugins/FFTJetProducer.h"
0029 
0030 // Framework include files
0031 #include "FWCore/Framework/interface/MakerMacros.h"
0032 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0033 
0034 // Data formats
0035 #include "DataFormats/Common/interface/View.h"
0036 #include "DataFormats/Common/interface/Handle.h"
0037 #include "DataFormats/VertexReco/interface/Vertex.h"
0038 #include "DataFormats/Candidate/interface/Candidate.h"
0039 
0040 #include "RecoJets/JetProducers/interface/JetSpecific.h"
0041 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0042 #include "Geometry/Records/interface/HcalRecNumberingRecord.h"
0043 
0044 // Loader for the lookup tables
0045 #include "JetMETCorrections/FFTJetObjects/interface/FFTJetLookupTableSequenceLoader.h"
0046 
0047 #define make_param(type, varname) const type& varname(ps.getParameter<type>(#varname))
0048 
0049 #define init_param(type, varname) varname(ps.getParameter<type>(#varname))
0050 
0051 // A generic switch statement based on jet type.
0052 // Defining it in a single place avoids potential errors
0053 // in case another jet type is introduced in the future.
0054 //
0055 // JPTJet is omitted for now: there is no reco::writeSpecific method
0056 // for it (see header JetSpecific.h in the JetProducers package)
0057 //
0058 #define jet_type_switch(method, arg1, arg2)                    \
0059   do {                                                         \
0060     switch (jetType) {                                         \
0061       case CALOJET:                                            \
0062         method<reco::CaloJet>(arg1, arg2);                     \
0063         break;                                                 \
0064       case PFJET:                                              \
0065         method<reco::PFJet>(arg1, arg2);                       \
0066         break;                                                 \
0067       case GENJET:                                             \
0068         method<reco::GenJet>(arg1, arg2);                      \
0069         break;                                                 \
0070       case TRACKJET:                                           \
0071         method<reco::TrackJet>(arg1, arg2);                    \
0072         break;                                                 \
0073       case BASICJET:                                           \
0074         method<reco::BasicJet>(arg1, arg2);                    \
0075         break;                                                 \
0076       default:                                                 \
0077         assert(!"ERROR in FFTJetProducer : invalid jet type."\
0078                " This is a bug. Please report."); \
0079     }                                                          \
0080   } while (0);
0081 
0082 namespace {
0083   struct LocalSortByPt {
0084     template <class Jet>
0085     inline bool operator()(const Jet& l, const Jet& r) const {
0086       return l.pt() > r.pt();
0087     }
0088   };
0089 }  // namespace
0090 
0091 using namespace fftjetcms;
0092 
0093 FFTJetProducer::Resolution FFTJetProducer::parse_resolution(const std::string& name) {
0094   if (!name.compare("fixed"))
0095     return FIXED;
0096   else if (!name.compare("maximallyStable"))
0097     return MAXIMALLY_STABLE;
0098   else if (!name.compare("globallyAdaptive"))
0099     return GLOBALLY_ADAPTIVE;
0100   else if (!name.compare("locallyAdaptive"))
0101     return LOCALLY_ADAPTIVE;
0102   else if (!name.compare("fromGenJets"))
0103     return FROM_GENJETS;
0104   else
0105     throw cms::Exception("FFTJetBadConfig") << "Invalid resolution specification \"" << name << "\"" << std::endl;
0106 }
0107 
0108 template <typename T>
0109 void FFTJetProducer::makeProduces(const std::string& alias, const std::string& tag) {
0110   produces<std::vector<reco::FFTAnyJet<T> > >(tag).setBranchAlias(alias);
0111 }
0112 
0113 //
0114 // constructors and destructor
0115 //
0116 FFTJetProducer::FFTJetProducer(const edm::ParameterSet& ps)
0117     : FFTJetInterface(ps),
0118       myConfiguration(ps),
0119       init_param(edm::InputTag, treeLabel),
0120       init_param(bool, useGriddedAlgorithm),
0121       init_param(bool, reuseExistingGrid),
0122       init_param(unsigned, maxIterations),
0123       init_param(unsigned, nJetsRequiredToConverge),
0124       init_param(double, convergenceDistance),
0125       init_param(bool, assignConstituents),
0126       init_param(bool, resumConstituents),
0127       init_param(bool, calculatePileup),
0128       init_param(bool, subtractPileup),
0129       init_param(bool, subtractPileupAs4Vec),
0130       init_param(edm::InputTag, pileupLabel),
0131       init_param(double, fixedScale),
0132       init_param(double, minStableScale),
0133       init_param(double, maxStableScale),
0134       init_param(double, stabilityAlpha),
0135       init_param(double, noiseLevel),
0136       init_param(unsigned, nClustersRequested),
0137       init_param(double, gridScanMaxEta),
0138       init_param(std::string, recombinationAlgorithm),
0139       init_param(bool, isCrisp),
0140       init_param(double, unlikelyBgWeight),
0141       init_param(double, recombinationDataCutoff),
0142       init_param(edm::InputTag, genJetsLabel),
0143       init_param(unsigned, maxInitialPreclusters),
0144       resolution(parse_resolution(ps.getParameter<std::string>("resolution"))),
0145       init_param(std::string, pileupTableRecord),
0146       init_param(std::string, pileupTableName),
0147       init_param(std::string, pileupTableCategory),
0148       init_param(bool, loadPileupFromDB),
0149 
0150       minLevel(0),
0151       maxLevel(0),
0152       usedLevel(0),
0153       unused(0.0),
0154       iterationsPerformed(1U),
0155       constituents(200) {
0156   // Check that the settings make sense
0157   if (resumConstituents && !assignConstituents)
0158     throw cms::Exception("FFTJetBadConfig") << "Can't resum constituents if they are not assigned" << std::endl;
0159 
0160   produces<reco::FFTJetProducerSummary>(outputLabel);
0161   const std::string alias(ps.getUntrackedParameter<std::string>("alias", outputLabel));
0162   jet_type_switch(makeProduces, alias, outputLabel);
0163 
0164   if (jetType == CALOJET) {
0165     geometry_token_ = esConsumes();
0166     topology_token_ = esConsumes();
0167   }
0168 
0169   // Build the set of pattern recognition scales.
0170   // This is needed in order to read the clustering tree
0171   // from the event record.
0172   iniScales = fftjet_ScaleSet_parser(ps.getParameter<edm::ParameterSet>("InitialScales"));
0173   checkConfig(iniScales, "invalid set of scales");
0174   std::sort(iniScales->begin(), iniScales->end(), std::greater<double>());
0175 
0176   if (storeInSinglePrecision())
0177     input_recotree_token_f_ = consumes<reco::PattRecoTree<float, reco::PattRecoPeak<float> > >(treeLabel);
0178   else
0179     input_recotree_token_d_ = consumes<reco::PattRecoTree<double, reco::PattRecoPeak<double> > >(treeLabel);
0180   input_genjet_token_ = consumes<std::vector<reco::FFTAnyJet<reco::GenJet> > >(genJetsLabel);
0181   input_energyflow_token_ = consumes<reco::DiscretizedEnergyFlow>(treeLabel);
0182   input_pusummary_token_ = consumes<reco::FFTJetPileupSummary>(pileupLabel);
0183 
0184   // Most of the configuration has to be performed inside
0185   // the "beginStream" method. This is because chaining of the
0186   // parsers between this base class and the derived classes
0187   // can not work from the constructor of the base class.
0188 }
0189 
0190 FFTJetProducer::~FFTJetProducer() {}
0191 
0192 //
0193 // member functions
0194 //
0195 template <class Real>
0196 void FFTJetProducer::loadSparseTreeData(
0197     const edm::Event& iEvent, const edm::EDGetTokenT<reco::PattRecoTree<Real, reco::PattRecoPeak<Real> > >& tok) {
0198   typedef reco::PattRecoTree<Real, reco::PattRecoPeak<Real> > StoredTree;
0199 
0200   // Get the input
0201   edm::Handle<StoredTree> input;
0202   iEvent.getByToken(tok, input);
0203 
0204   if (!input->isSparse())
0205     throw cms::Exception("FFTJetBadConfig") << "The stored clustering tree is not sparse" << std::endl;
0206 
0207   sparsePeakTreeFromStorable(*input, iniScales.get(), getEventScale(), &sparseTree);
0208   sparseTree.sortNodes();
0209   fftjet::updateSplitMergeTimes(sparseTree, sparseTree.minScale(), sparseTree.maxScale());
0210 }
0211 
0212 void FFTJetProducer::genJetPreclusters(const SparseTree& /* tree */,
0213                                        edm::Event& iEvent,
0214                                        const edm::EventSetup& /* iSetup */,
0215                                        const fftjet::Functor1<bool, fftjet::Peak>& peakSelect,
0216                                        std::vector<fftjet::Peak>* preclusters) {
0217   typedef reco::FFTAnyJet<reco::GenJet> InputJet;
0218   typedef std::vector<InputJet> InputCollection;
0219 
0220   edm::Handle<InputCollection> input;
0221   iEvent.getByToken(input_genjet_token_, input);
0222 
0223   const unsigned sz = input->size();
0224   preclusters->reserve(sz);
0225   for (unsigned i = 0; i < sz; ++i) {
0226     const RecoFFTJet& jet(jetFromStorable((*input)[i].getFFTSpecific()));
0227     fftjet::Peak p(jet.precluster());
0228     const double scale(p.scale());
0229     p.setEtaPhi(jet.vec().Eta(), jet.vec().Phi());
0230     p.setMagnitude(jet.vec().Pt() / scale / scale);
0231     p.setStatus(resolution);
0232     if (peakSelect(p))
0233       preclusters->push_back(p);
0234   }
0235 }
0236 
0237 void FFTJetProducer::selectPreclusters(const SparseTree& tree,
0238                                        const fftjet::Functor1<bool, fftjet::Peak>& peakSelect,
0239                                        std::vector<fftjet::Peak>* preclusters) {
0240   nodes.clear();
0241   selectTreeNodes(tree, peakSelect, &nodes);
0242 
0243   // Fill out the vector of preclusters using the tree node ids
0244   const unsigned nNodes = nodes.size();
0245   const SparseTree::NodeId* pnodes = nNodes ? &nodes[0] : nullptr;
0246   preclusters->reserve(nNodes);
0247   for (unsigned i = 0; i < nNodes; ++i)
0248     preclusters->push_back(sparseTree.uncheckedNode(pnodes[i]).getCluster());
0249 
0250   // Remember the node id in the precluster and set
0251   // the status word to indicate the resolution scheme used
0252   fftjet::Peak* clusters = nNodes ? &(*preclusters)[0] : nullptr;
0253   for (unsigned i = 0; i < nNodes; ++i) {
0254     clusters[i].setCode(pnodes[i]);
0255     clusters[i].setStatus(resolution);
0256   }
0257 }
0258 
0259 void FFTJetProducer::selectTreeNodes(const SparseTree& tree,
0260                                      const fftjet::Functor1<bool, fftjet::Peak>& peakSelect,
0261                                      std::vector<SparseTree::NodeId>* mynodes) {
0262   minLevel = maxLevel = usedLevel = 0;
0263 
0264   // Get the tree nodes which pass the cuts
0265   // (according to the selected resolution strategy)
0266   switch (resolution) {
0267     case FIXED: {
0268       usedLevel = tree.getLevel(fixedScale);
0269       tree.getPassingNodes(usedLevel, peakSelect, mynodes);
0270     } break;
0271 
0272     case MAXIMALLY_STABLE: {
0273       const unsigned minStartingLevel = maxStableScale > 0.0 ? tree.getLevel(maxStableScale) : 0;
0274       const unsigned maxStartingLevel = minStableScale > 0.0 ? tree.getLevel(minStableScale) : UINT_MAX;
0275 
0276       if (tree.stableClusterCount(
0277               peakSelect, &minLevel, &maxLevel, stabilityAlpha, minStartingLevel, maxStartingLevel)) {
0278         usedLevel = (minLevel + maxLevel) / 2;
0279         tree.getPassingNodes(usedLevel, peakSelect, mynodes);
0280       }
0281     } break;
0282 
0283     case GLOBALLY_ADAPTIVE: {
0284       const bool stable = tree.clusterCountLevels(nClustersRequested, peakSelect, &minLevel, &maxLevel);
0285       if (minLevel || maxLevel) {
0286         usedLevel = (minLevel + maxLevel) / 2;
0287         if (!stable) {
0288           const int maxlev = tree.maxStoredLevel();
0289           bool levelFound = false;
0290           for (int delta = 0; delta <= maxlev && !levelFound; ++delta)
0291             for (int ifac = 1; ifac > -2 && !levelFound; ifac -= 2) {
0292               const int level = usedLevel + ifac * delta;
0293               if (level > 0 && level <= maxlev)
0294                 if (occupancy[level] == nClustersRequested) {
0295                   usedLevel = level;
0296                   levelFound = true;
0297                 }
0298             }
0299           assert(levelFound);
0300         }
0301       } else {
0302         // Can't find that exact number of preclusters.
0303         // Try to get the next best thing.
0304         usedLevel = 1;
0305         const unsigned occ1 = occupancy[1];
0306         if (nClustersRequested >= occ1) {
0307           const unsigned maxlev = tree.maxStoredLevel();
0308           if (nClustersRequested > occupancy[maxlev])
0309             usedLevel = maxlev;
0310           else {
0311             // It would be nice to use "lower_bound" here,
0312             // but the occupancy is not necessarily monotonous.
0313             unsigned bestDelta = nClustersRequested > occ1 ? nClustersRequested - occ1 : occ1 - nClustersRequested;
0314             for (unsigned level = 2; level <= maxlev; ++level) {
0315               const unsigned n = occupancy[level];
0316               const unsigned d = nClustersRequested > n ? nClustersRequested - n : n - nClustersRequested;
0317               if (d < bestDelta) {
0318                 bestDelta = d;
0319                 usedLevel = level;
0320               }
0321             }
0322           }
0323         }
0324       }
0325       tree.getPassingNodes(usedLevel, peakSelect, mynodes);
0326     } break;
0327 
0328     case LOCALLY_ADAPTIVE: {
0329       usedLevel = tree.getLevel(fixedScale);
0330       tree.getMagS2OptimalNodes(peakSelect, nClustersRequested, usedLevel, mynodes, &thresholds);
0331     } break;
0332 
0333     default:
0334       assert(!"ERROR in FFTJetProducer::selectTreeNodes : "
0335                "should never get here! This is a bug. Please report.");
0336   }
0337 }
0338 
0339 void FFTJetProducer::prepareRecombinationScales() {
0340   const unsigned nClus = preclusters.size();
0341   if (nClus) {
0342     fftjet::Peak* clus = &preclusters[0];
0343     fftjet::Functor1<double, fftjet::Peak>& scaleCalc(*recoScaleCalcPeak);
0344     fftjet::Functor1<double, fftjet::Peak>& ratioCalc(*recoScaleRatioCalcPeak);
0345     fftjet::Functor1<double, fftjet::Peak>& factorCalc(*memberFactorCalcPeak);
0346 
0347     for (unsigned i = 0; i < nClus; ++i) {
0348       clus[i].setRecoScale(scaleCalc(clus[i]));
0349       clus[i].setRecoScaleRatio(ratioCalc(clus[i]));
0350       clus[i].setMembershipFactor(factorCalc(clus[i]));
0351     }
0352   }
0353 }
0354 
0355 void FFTJetProducer::buildGridAlg() {
0356   int minBin = energyFlow->getEtaBin(-gridScanMaxEta);
0357   if (minBin < 0)
0358     minBin = 0;
0359   int maxBin = energyFlow->getEtaBin(gridScanMaxEta) + 1;
0360   if (maxBin < 0)
0361     maxBin = 0;
0362 
0363   fftjet::DefaultRecombinationAlgFactory<Real, VectorLike, BgData, VBuilder> factory;
0364   if (factory[recombinationAlgorithm] == nullptr)
0365     throw cms::Exception("FFTJetBadConfig")
0366         << "Invalid grid recombination algorithm \"" << recombinationAlgorithm << "\"" << std::endl;
0367   gridAlg = std::unique_ptr<GridAlg>(factory[recombinationAlgorithm]->create(jetMembershipFunction.get(),
0368                                                                              bgMembershipFunction.get(),
0369                                                                              unlikelyBgWeight,
0370                                                                              recombinationDataCutoff,
0371                                                                              isCrisp,
0372                                                                              false,
0373                                                                              assignConstituents,
0374                                                                              minBin,
0375                                                                              maxBin));
0376 }
0377 
0378 bool FFTJetProducer::loadEnergyFlow(const edm::Event& iEvent, std::unique_ptr<fftjet::Grid2d<fftjetcms::Real> >& flow) {
0379   edm::Handle<reco::DiscretizedEnergyFlow> input;
0380   iEvent.getByToken(input_energyflow_token_, input);
0381 
0382   // Make sure that the grid is compatible with the stored one
0383   bool rebuildGrid = flow.get() == nullptr;
0384   if (!rebuildGrid)
0385     rebuildGrid =
0386         !(flow->nEta() == input->nEtaBins() && flow->nPhi() == input->nPhiBins() && flow->etaMin() == input->etaMin() &&
0387           flow->etaMax() == input->etaMax() && flow->phiBin0Edge() == input->phiBin0Edge());
0388   if (rebuildGrid) {
0389     // We should not get here very often...
0390     flow = std::make_unique<fftjet::Grid2d<Real> >(
0391         input->nEtaBins(), input->etaMin(), input->etaMax(), input->nPhiBins(), input->phiBin0Edge(), input->title());
0392   }
0393   flow->blockSet(input->data(), input->nEtaBins(), input->nPhiBins());
0394   return rebuildGrid;
0395 }
0396 
0397 bool FFTJetProducer::checkConvergence(const std::vector<RecoFFTJet>& previous, std::vector<RecoFFTJet>& nextSet) {
0398   fftjet::Functor2<double, RecoFFTJet, RecoFFTJet>& distanceCalc(*jetDistanceCalc);
0399 
0400   const unsigned nJets = previous.size();
0401   if (nJets != nextSet.size())
0402     return false;
0403 
0404   const RecoFFTJet* prev = &previous[0];
0405   RecoFFTJet* next = &nextSet[0];
0406 
0407   // Calculate convergence distances for all jets
0408   bool converged = true;
0409   for (unsigned i = 0; i < nJets; ++i) {
0410     const double d = distanceCalc(prev[i], next[i]);
0411     next[i].setConvergenceDistance(d);
0412     if (i < nJetsRequiredToConverge && d > convergenceDistance)
0413       converged = false;
0414   }
0415 
0416   return converged;
0417 }
0418 
0419 unsigned FFTJetProducer::iterateJetReconstruction() {
0420   fftjet::Functor1<double, RecoFFTJet>& scaleCalc(*recoScaleCalcJet);
0421   fftjet::Functor1<double, RecoFFTJet>& ratioCalc(*recoScaleRatioCalcJet);
0422   fftjet::Functor1<double, RecoFFTJet>& factorCalc(*memberFactorCalcJet);
0423 
0424   unsigned nJets = recoJets.size();
0425   unsigned iterNum = 1U;
0426   bool converged = false;
0427   for (; iterNum < maxIterations && !converged; ++iterNum) {
0428     // Recreate the vector of preclusters using the jets
0429     const RecoFFTJet* jets = &recoJets[0];
0430     iterPreclusters.clear();
0431     iterPreclusters.reserve(nJets);
0432     for (unsigned i = 0; i < nJets; ++i) {
0433       const RecoFFTJet& jet(jets[i]);
0434       fftjet::Peak p(jet.precluster());
0435       p.setEtaPhi(jet.vec().Eta(), jet.vec().Phi());
0436       p.setRecoScale(scaleCalc(jet));
0437       p.setRecoScaleRatio(ratioCalc(jet));
0438       p.setMembershipFactor(factorCalc(jet));
0439       iterPreclusters.push_back(p);
0440     }
0441 
0442     // Run the algorithm
0443     int status = 0;
0444     if (useGriddedAlgorithm)
0445       status = gridAlg->run(iterPreclusters, *energyFlow, &noiseLevel, 1U, 1U, &iterJets, &unclustered, &unused);
0446     else
0447       status = recoAlg->run(iterPreclusters, eventData, &noiseLevel, 1U, &iterJets, &unclustered, &unused);
0448     if (status)
0449       throw cms::Exception("FFTJetInterface") << "FFTJet algorithm failed" << std::endl;
0450 
0451     // As it turns out, it is possible, in very rare cases,
0452     // to have iterJets.size() != nJets at this point
0453 
0454     // Figure out if the iterations have converged
0455     converged = checkConvergence(recoJets, iterJets);
0456 
0457     // Prepare for the next cycle
0458     iterJets.swap(recoJets);
0459     nJets = recoJets.size();
0460   }
0461 
0462   // Check that we have the correct number of preclusters
0463   if (preclusters.size() != nJets) {
0464     assert(nJets < preclusters.size());
0465     removeFakePreclusters();
0466     assert(preclusters.size() == nJets);
0467   }
0468 
0469   // Plug in the original precluster coordinates into the result
0470   RecoFFTJet* jets = &recoJets[0];
0471   for (unsigned i = 0; i < nJets; ++i) {
0472     const fftjet::Peak& oldp(preclusters[i]);
0473     jets[i].setPeakEtaPhi(oldp.eta(), oldp.phi());
0474   }
0475 
0476   // If we have converged on the last cycle, the result
0477   // would be indistinguishable from no convergence.
0478   // Because of this, raise the counter by one to indicate
0479   // the case when the convergence is not achieved.
0480   if (!converged)
0481     ++iterNum;
0482 
0483   return iterNum;
0484 }
0485 
0486 void FFTJetProducer::determineGriddedConstituents() {
0487   const unsigned nJets = recoJets.size();
0488   const unsigned* clusterMask = gridAlg->getClusterMask();
0489   const int nEta = gridAlg->getLastNEta();
0490   const int nPhi = gridAlg->getLastNPhi();
0491   const fftjet::Grid2d<Real>& g(*energyFlow);
0492 
0493   const unsigned nInputs = eventData.size();
0494   const VectorLike* inp = nInputs ? &eventData[0] : nullptr;
0495   const unsigned* candIdx = nInputs ? &candidateIndex[0] : nullptr;
0496   for (unsigned i = 0; i < nInputs; ++i) {
0497     const VectorLike& item(inp[i]);
0498     const int iPhi = g.getPhiBin(item.Phi());
0499     const int iEta = g.getEtaBin(item.Eta());
0500     const unsigned mask = iEta >= 0 && iEta < nEta ? clusterMask[iEta * nPhi + iPhi] : 0;
0501     assert(mask <= nJets);
0502     constituents[mask].push_back(inputCollection->ptrAt(candIdx[i]));
0503   }
0504 }
0505 
0506 void FFTJetProducer::determineVectorConstituents() {
0507   const unsigned nJets = recoJets.size();
0508   const unsigned* clusterMask = recoAlg->getClusterMask();
0509   const unsigned maskLength = recoAlg->getLastNData();
0510   assert(maskLength == eventData.size());
0511 
0512   const unsigned* candIdx = maskLength ? &candidateIndex[0] : nullptr;
0513   for (unsigned i = 0; i < maskLength; ++i) {
0514     // In FFTJet, the mask value of 0 corresponds to unclustered
0515     // energy. We will do the same here. Jet numbers are therefore
0516     // shifted by 1 wrt constituents vector, and constituents[1]
0517     // corresponds to jet number 0.
0518     const unsigned mask = clusterMask[i];
0519     assert(mask <= nJets);
0520     constituents[mask].push_back(inputCollection->ptrAt(candIdx[i]));
0521   }
0522 }
0523 
0524 // The following code more or less coincides with the similar method
0525 // implemented in VirtualJetProducer
0526 template <typename T>
0527 void FFTJetProducer::writeJets(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0528   using namespace reco;
0529 
0530   typedef FFTAnyJet<T> OutputJet;
0531   typedef std::vector<OutputJet> OutputCollection;
0532 
0533   // Area of a single eta-phi cell for jet area calculations.
0534   // Set it to 0 in case the module configuration does not allow
0535   // us to calculate jet areas reliably.
0536   double cellArea = useGriddedAlgorithm && recombinationDataCutoff < 0.0
0537                         ? energyFlow->etaBinWidth() * energyFlow->phiBinWidth()
0538                         : 0.0;
0539 
0540   if (calculatePileup)
0541     cellArea = pileupEnergyFlow->etaBinWidth() * pileupEnergyFlow->phiBinWidth();
0542 
0543   // allocate output jet collection
0544   auto jets = std::make_unique<OutputCollection>();
0545   const unsigned nJets = recoJets.size();
0546   jets->reserve(nJets);
0547 
0548   bool sorted = true;
0549   double previousPt = DBL_MAX;
0550   for (unsigned ijet = 0; ijet < nJets; ++ijet) {
0551     RecoFFTJet& myjet(recoJets[ijet]);
0552 
0553     // Check if we should resum jet constituents
0554     VectorLike jet4vec(myjet.vec());
0555     if (resumConstituents) {
0556       VectorLike sum(0.0, 0.0, 0.0, 0.0);
0557       const unsigned nCon = constituents[ijet + 1].size();
0558       const reco::CandidatePtr* cn = nCon ? &constituents[ijet + 1][0] : nullptr;
0559       for (unsigned i = 0; i < nCon; ++i)
0560         sum += cn[i]->p4();
0561       jet4vec = sum;
0562       setJetStatusBit(&myjet, CONSTITUENTS_RESUMMED, true);
0563     }
0564 
0565     // Subtract the pile-up
0566     if (calculatePileup && subtractPileup) {
0567       jet4vec = adjustForPileup(jet4vec, pileup[ijet], subtractPileupAs4Vec);
0568       if (subtractPileupAs4Vec)
0569         setJetStatusBit(&myjet, PILEUP_SUBTRACTED_4VEC, true);
0570       else
0571         setJetStatusBit(&myjet, PILEUP_SUBTRACTED_PT, true);
0572     }
0573 
0574     // Write the specifics to the jet (simultaneously sets 4-vector,
0575     // vertex, constituents). These are overridden functions that will
0576     // call the appropriate specific code.
0577     T jet;
0578     if constexpr (std::is_same_v<T, reco::CaloJet>) {
0579       const CaloGeometry& geometry = iSetup.getData(geometry_token_);
0580       const HcalTopology& topology = iSetup.getData(topology_token_);
0581       writeSpecific(jet, jet4vec, vertexUsed(), constituents[ijet + 1], geometry, topology);
0582     } else {
0583       writeSpecific(jet, jet4vec, vertexUsed(), constituents[ijet + 1]);
0584     }
0585 
0586     // calcuate the jet area
0587     double ncells = myjet.ncells();
0588     if (calculatePileup) {
0589       ncells = cellCountsVec[ijet];
0590       setJetStatusBit(&myjet, PILEUP_CALCULATED, true);
0591     }
0592     jet.setJetArea(cellArea * ncells);
0593 
0594     // add jet to the list
0595     FFTJet<float> fj(jetToStorable<float>(myjet));
0596     fj.setFourVec(jet4vec);
0597     if (calculatePileup) {
0598       fj.setPileup(pileup[ijet]);
0599       fj.setNCells(ncells);
0600     }
0601     jets->push_back(OutputJet(jet, fj));
0602 
0603     // Check whether the sequence remains sorted by pt
0604     const double pt = jet.pt();
0605     if (pt > previousPt)
0606       sorted = false;
0607     previousPt = pt;
0608   }
0609 
0610   // Sort the collection
0611   if (!sorted)
0612     std::sort(jets->begin(), jets->end(), LocalSortByPt());
0613 
0614   // put the collection into the event
0615   iEvent.put(std::move(jets), outputLabel);
0616 }
0617 
0618 void FFTJetProducer::saveResults(edm::Event& ev, const edm::EventSetup& iSetup, const unsigned nPreclustersFound) {
0619   // Write recombined jets
0620   jet_type_switch(writeJets, ev, iSetup);
0621 
0622   // Check if we should resum unclustered energy constituents
0623   VectorLike unclusE(unclustered);
0624   if (resumConstituents) {
0625     VectorLike sum(0.0, 0.0, 0.0, 0.0);
0626     const unsigned nCon = constituents[0].size();
0627     const reco::CandidatePtr* cn = nCon ? &constituents[0][0] : nullptr;
0628     for (unsigned i = 0; i < nCon; ++i)
0629       sum += cn[i]->p4();
0630     unclusE = sum;
0631   }
0632 
0633   // Write the jet reconstruction summary
0634   const double minScale = minLevel ? sparseTree.getScale(minLevel) : 0.0;
0635   const double maxScale = maxLevel ? sparseTree.getScale(maxLevel) : 0.0;
0636   const double scaleUsed = usedLevel ? sparseTree.getScale(usedLevel) : 0.0;
0637 
0638   ev.put(
0639       std::make_unique<reco::FFTJetProducerSummary>(thresholds,
0640                                                     occupancy,
0641                                                     unclusE,
0642                                                     constituents[0],
0643                                                     unused,
0644                                                     minScale,
0645                                                     maxScale,
0646                                                     scaleUsed,
0647                                                     nPreclustersFound,
0648                                                     iterationsPerformed,
0649                                                     iterationsPerformed == 1U || iterationsPerformed <= maxIterations),
0650       outputLabel);
0651 }
0652 
0653 // ------------ method called to for each event  ------------
0654 void FFTJetProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0655   // Load the clustering tree made by FFTJetPatRecoProducer
0656   if (storeInSinglePrecision())
0657     loadSparseTreeData<float>(iEvent, input_recotree_token_f_);
0658   else
0659     loadSparseTreeData<double>(iEvent, input_recotree_token_d_);
0660 
0661   // Do we need to load the candidate collection?
0662   if (assignConstituents || !(useGriddedAlgorithm && reuseExistingGrid))
0663     loadInputCollection(iEvent);
0664 
0665   // Do we need to have discretized energy flow?
0666   if (useGriddedAlgorithm) {
0667     if (reuseExistingGrid) {
0668       if (loadEnergyFlow(iEvent, energyFlow))
0669         buildGridAlg();
0670     } else
0671       discretizeEnergyFlow();
0672   }
0673 
0674   // Calculate cluster occupancy as a function of level number
0675   sparseTree.occupancyInScaleSpace(*peakSelector, &occupancy);
0676 
0677   // Select the preclusters using the requested resolution scheme
0678   preclusters.clear();
0679   if (resolution == FROM_GENJETS)
0680     genJetPreclusters(sparseTree, iEvent, iSetup, *peakSelector, &preclusters);
0681   else
0682     selectPreclusters(sparseTree, *peakSelector, &preclusters);
0683   if (preclusters.size() > maxInitialPreclusters) {
0684     std::sort(preclusters.begin(), preclusters.end(), std::greater<fftjet::Peak>());
0685     preclusters.erase(preclusters.begin() + maxInitialPreclusters, preclusters.end());
0686   }
0687 
0688   // Prepare to run the jet recombination procedure
0689   prepareRecombinationScales();
0690 
0691   // Assign membership functions to preclusters. If this function
0692   // is not overriden in a derived class, default algorithm membership
0693   // function will be used for every cluster.
0694   assignMembershipFunctions(&preclusters);
0695 
0696   // Count the preclusters going in
0697   unsigned nPreclustersFound = 0U;
0698   const unsigned npre = preclusters.size();
0699   for (unsigned i = 0; i < npre; ++i)
0700     if (preclusters[i].membershipFactor() > 0.0)
0701       ++nPreclustersFound;
0702 
0703   // Run the recombination algorithm once
0704   int status = 0;
0705   if (useGriddedAlgorithm)
0706     status = gridAlg->run(preclusters, *energyFlow, &noiseLevel, 1U, 1U, &recoJets, &unclustered, &unused);
0707   else
0708     status = recoAlg->run(preclusters, eventData, &noiseLevel, 1U, &recoJets, &unclustered, &unused);
0709   if (status)
0710     throw cms::Exception("FFTJetInterface") << "FFTJet algorithm failed (first iteration)" << std::endl;
0711 
0712   // If requested, iterate the jet recombination procedure
0713   if (maxIterations > 1U && !recoJets.empty()) {
0714     // It is possible to have a smaller number of jets than we had
0715     // preclusters. Fake preclusters are possible, but for a good
0716     // choice of pattern recognition kernel their presence should
0717     // be infrequent. However, any fake preclusters will throw the
0718     // iterative reconstruction off balance. Deal with the problem now.
0719     const unsigned nJets = recoJets.size();
0720     if (preclusters.size() != nJets) {
0721       assert(nJets < preclusters.size());
0722       removeFakePreclusters();
0723     }
0724     iterationsPerformed = iterateJetReconstruction();
0725   } else
0726     iterationsPerformed = 1U;
0727 
0728   // Determine jet constituents. FFTJet returns a map
0729   // of constituents which is inverse to what we need here.
0730   const unsigned nJets = recoJets.size();
0731   if (constituents.size() <= nJets)
0732     constituents.resize(nJets + 1U);
0733   if (assignConstituents) {
0734     for (unsigned i = 0; i <= nJets; ++i)
0735       constituents[i].clear();
0736     if (useGriddedAlgorithm)
0737       determineGriddedConstituents();
0738     else
0739       determineVectorConstituents();
0740   }
0741 
0742   // Figure out the pile-up
0743   if (calculatePileup) {
0744     if (loadPileupFromDB)
0745       determinePileupDensityFromDB(iEvent, iSetup, pileupEnergyFlow);
0746     else
0747       determinePileupDensityFromConfig(iEvent, pileupEnergyFlow);
0748     determinePileup();
0749     assert(pileup.size() == recoJets.size());
0750   }
0751 
0752   // Write out the results
0753   saveResults(iEvent, iSetup, nPreclustersFound);
0754 }
0755 
0756 std::unique_ptr<fftjet::Functor1<bool, fftjet::Peak> > FFTJetProducer::parse_peakSelector(const edm::ParameterSet& ps) {
0757   return fftjet_PeakSelector_parser(ps.getParameter<edm::ParameterSet>("PeakSelectorConfiguration"));
0758 }
0759 
0760 // Parser for the jet membership function
0761 std::unique_ptr<fftjet::ScaleSpaceKernel> FFTJetProducer::parse_jetMembershipFunction(const edm::ParameterSet& ps) {
0762   return fftjet_MembershipFunction_parser(ps.getParameter<edm::ParameterSet>("jetMembershipFunction"));
0763 }
0764 
0765 // Parser for the background membership function
0766 std::unique_ptr<AbsBgFunctor> FFTJetProducer::parse_bgMembershipFunction(const edm::ParameterSet& ps) {
0767   return fftjet_BgFunctor_parser(ps.getParameter<edm::ParameterSet>("bgMembershipFunction"));
0768 }
0769 
0770 // Calculator for the recombination scale
0771 std::unique_ptr<fftjet::Functor1<double, fftjet::Peak> > FFTJetProducer::parse_recoScaleCalcPeak(
0772     const edm::ParameterSet& ps) {
0773   return fftjet_PeakFunctor_parser(ps.getParameter<edm::ParameterSet>("recoScaleCalcPeak"));
0774 }
0775 
0776 // Pile-up density calculator
0777 std::unique_ptr<fftjetcms::AbsPileupCalculator> FFTJetProducer::parse_pileupDensityCalc(const edm::ParameterSet& ps) {
0778   return fftjet_PileupCalculator_parser(ps.getParameter<edm::ParameterSet>("pileupDensityCalc"));
0779 }
0780 
0781 // Calculator for the recombination scale ratio
0782 std::unique_ptr<fftjet::Functor1<double, fftjet::Peak> > FFTJetProducer::parse_recoScaleRatioCalcPeak(
0783     const edm::ParameterSet& ps) {
0784   return fftjet_PeakFunctor_parser(ps.getParameter<edm::ParameterSet>("recoScaleRatioCalcPeak"));
0785 }
0786 
0787 // Calculator for the membership function factor
0788 std::unique_ptr<fftjet::Functor1<double, fftjet::Peak> > FFTJetProducer::parse_memberFactorCalcPeak(
0789     const edm::ParameterSet& ps) {
0790   return fftjet_PeakFunctor_parser(ps.getParameter<edm::ParameterSet>("memberFactorCalcPeak"));
0791 }
0792 
0793 std::unique_ptr<fftjet::Functor1<double, FFTJetProducer::RecoFFTJet> > FFTJetProducer::parse_recoScaleCalcJet(
0794     const edm::ParameterSet& ps) {
0795   return fftjet_JetFunctor_parser(ps.getParameter<edm::ParameterSet>("recoScaleCalcJet"));
0796 }
0797 
0798 std::unique_ptr<fftjet::Functor1<double, FFTJetProducer::RecoFFTJet> > FFTJetProducer::parse_recoScaleRatioCalcJet(
0799     const edm::ParameterSet& ps) {
0800   return fftjet_JetFunctor_parser(ps.getParameter<edm::ParameterSet>("recoScaleRatioCalcJet"));
0801 }
0802 
0803 std::unique_ptr<fftjet::Functor1<double, FFTJetProducer::RecoFFTJet> > FFTJetProducer::parse_memberFactorCalcJet(
0804     const edm::ParameterSet& ps) {
0805   return fftjet_JetFunctor_parser(ps.getParameter<edm::ParameterSet>("memberFactorCalcJet"));
0806 }
0807 
0808 std::unique_ptr<fftjet::Functor2<double, FFTJetProducer::RecoFFTJet, FFTJetProducer::RecoFFTJet> >
0809 FFTJetProducer::parse_jetDistanceCalc(const edm::ParameterSet& ps) {
0810   return fftjet_JetDistance_parser(ps.getParameter<edm::ParameterSet>("jetDistanceCalc"));
0811 }
0812 
0813 void FFTJetProducer::assignMembershipFunctions(std::vector<fftjet::Peak>*) {}
0814 
0815 // ------------ method called once each job just before starting event loop
0816 void FFTJetProducer::beginStream(edm::StreamID) {
0817   const edm::ParameterSet& ps(myConfiguration);
0818 
0819   // Parse the peak selector definition
0820   peakSelector = parse_peakSelector(ps);
0821   checkConfig(peakSelector, "invalid peak selector");
0822 
0823   jetMembershipFunction = parse_jetMembershipFunction(ps);
0824   checkConfig(jetMembershipFunction, "invalid jet membership function");
0825 
0826   bgMembershipFunction = parse_bgMembershipFunction(ps);
0827   checkConfig(bgMembershipFunction, "invalid noise membership function");
0828 
0829   // Build the energy recombination algorithm
0830   if (!useGriddedAlgorithm) {
0831     fftjet::DefaultVectorRecombinationAlgFactory<VectorLike, BgData, VBuilder> factory;
0832     if (factory[recombinationAlgorithm] == nullptr)
0833       throw cms::Exception("FFTJetBadConfig")
0834           << "Invalid vector recombination algorithm \"" << recombinationAlgorithm << "\"" << std::endl;
0835     recoAlg = std::unique_ptr<RecoAlg>(factory[recombinationAlgorithm]->create(jetMembershipFunction.get(),
0836                                                                                &VectorLike::Et,
0837                                                                                &VectorLike::Eta,
0838                                                                                &VectorLike::Phi,
0839                                                                                bgMembershipFunction.get(),
0840                                                                                unlikelyBgWeight,
0841                                                                                isCrisp,
0842                                                                                false,
0843                                                                                assignConstituents));
0844   } else if (!reuseExistingGrid) {
0845     energyFlow = fftjet_Grid2d_parser(ps.getParameter<edm::ParameterSet>("GridConfiguration"));
0846     checkConfig(energyFlow, "invalid discretization grid");
0847     buildGridAlg();
0848   }
0849 
0850   // Create the grid subsequently used for pile-up subtraction
0851   if (calculatePileup) {
0852     pileupEnergyFlow = fftjet_Grid2d_parser(ps.getParameter<edm::ParameterSet>("PileupGridConfiguration"));
0853     checkConfig(pileupEnergyFlow, "invalid pileup density grid");
0854 
0855     if (!loadPileupFromDB) {
0856       pileupDensityCalc = parse_pileupDensityCalc(ps);
0857       checkConfig(pileupDensityCalc, "invalid pile-up density calculator");
0858     }
0859   }
0860 
0861   // Parse the calculator of the recombination scale
0862   recoScaleCalcPeak = parse_recoScaleCalcPeak(ps);
0863   checkConfig(recoScaleCalcPeak,
0864               "invalid spec for the "
0865               "reconstruction scale calculator from peaks");
0866 
0867   // Parse the calculator of the recombination scale ratio
0868   recoScaleRatioCalcPeak = parse_recoScaleRatioCalcPeak(ps);
0869   checkConfig(recoScaleRatioCalcPeak,
0870               "invalid spec for the "
0871               "reconstruction scale ratio calculator from peaks");
0872 
0873   // Calculator for the membership function factor
0874   memberFactorCalcPeak = parse_memberFactorCalcPeak(ps);
0875   checkConfig(memberFactorCalcPeak,
0876               "invalid spec for the "
0877               "membership function factor calculator from peaks");
0878 
0879   if (maxIterations > 1) {
0880     // We are going to run iteratively. Make required objects.
0881     recoScaleCalcJet = parse_recoScaleCalcJet(ps);
0882     checkConfig(recoScaleCalcJet,
0883                 "invalid spec for the "
0884                 "reconstruction scale calculator from jets");
0885 
0886     recoScaleRatioCalcJet = parse_recoScaleRatioCalcJet(ps);
0887     checkConfig(recoScaleRatioCalcJet,
0888                 "invalid spec for the "
0889                 "reconstruction scale ratio calculator from jets");
0890 
0891     memberFactorCalcJet = parse_memberFactorCalcJet(ps);
0892     checkConfig(memberFactorCalcJet,
0893                 "invalid spec for the "
0894                 "membership function factor calculator from jets");
0895 
0896     jetDistanceCalc = parse_jetDistanceCalc(ps);
0897     checkConfig(memberFactorCalcJet,
0898                 "invalid spec for the "
0899                 "jet distance calculator");
0900   }
0901 }
0902 
0903 void FFTJetProducer::removeFakePreclusters() {
0904   // There are two possible reasons for fake preclusters:
0905   // 1. Membership factor was set to 0
0906   // 2. Genuine problem with pattern recognition
0907   //
0908   // Anyway, we need to match jets to preclusters and keep
0909   // only those preclusters that have been matched
0910   //
0911   std::vector<int> matchTable;
0912   const unsigned nmatched = matchOneToOne(recoJets, preclusters, JetToPeakDistance(), &matchTable);
0913 
0914   // Ensure that all jets have been matched.
0915   // If not, we must have a bug somewhere.
0916   assert(nmatched == recoJets.size());
0917 
0918   // Collect all matched preclusters
0919   iterPreclusters.clear();
0920   iterPreclusters.reserve(nmatched);
0921   for (unsigned i = 0; i < nmatched; ++i)
0922     iterPreclusters.push_back(preclusters[matchTable[i]]);
0923   iterPreclusters.swap(preclusters);
0924 }
0925 
0926 void FFTJetProducer::setJetStatusBit(RecoFFTJet* jet, const int mask, const bool value) {
0927   int status = jet->status();
0928   if (value)
0929     status |= mask;
0930   else
0931     status &= ~mask;
0932   jet->setStatus(status);
0933 }
0934 
0935 void FFTJetProducer::determinePileupDensityFromConfig(const edm::Event& iEvent,
0936                                                       std::unique_ptr<fftjet::Grid2d<fftjetcms::Real> >& density) {
0937   edm::Handle<reco::FFTJetPileupSummary> summary;
0938   iEvent.getByToken(input_pusummary_token_, summary);
0939 
0940   const reco::FFTJetPileupSummary& s(*summary);
0941   const AbsPileupCalculator& calc(*pileupDensityCalc);
0942   const bool phiDependent = calc.isPhiDependent();
0943 
0944   fftjet::Grid2d<Real>& g(*density);
0945   const unsigned nEta = g.nEta();
0946   const unsigned nPhi = g.nPhi();
0947 
0948   for (unsigned ieta = 0; ieta < nEta; ++ieta) {
0949     const double eta(g.etaBinCenter(ieta));
0950 
0951     if (phiDependent) {
0952       for (unsigned iphi = 0; iphi < nPhi; ++iphi) {
0953         const double phi(g.phiBinCenter(iphi));
0954         g.uncheckedSetBin(ieta, iphi, calc(eta, phi, s));
0955       }
0956     } else {
0957       const double pil = calc(eta, 0.0, s);
0958       for (unsigned iphi = 0; iphi < nPhi; ++iphi)
0959         g.uncheckedSetBin(ieta, iphi, pil);
0960     }
0961   }
0962 }
0963 
0964 void FFTJetProducer::determinePileupDensityFromDB(const edm::Event& iEvent,
0965                                                   const edm::EventSetup& iSetup,
0966                                                   std::unique_ptr<fftjet::Grid2d<fftjetcms::Real> >& density) {
0967   edm::ESHandle<FFTJetLookupTableSequence> h;
0968   StaticFFTJetLookupTableSequenceLoader::instance().load(iSetup, pileupTableRecord, h);
0969   std::shared_ptr<npstat::StorableMultivariateFunctor> f = (*h)[pileupTableCategory][pileupTableName];
0970 
0971   edm::Handle<reco::FFTJetPileupSummary> summary;
0972   iEvent.getByToken(input_pusummary_token_, summary);
0973 
0974   const float rho = summary->pileupRho();
0975   const bool phiDependent = f->minDim() == 3U;
0976 
0977   fftjet::Grid2d<Real>& g(*density);
0978   const unsigned nEta = g.nEta();
0979   const unsigned nPhi = g.nPhi();
0980 
0981   double functorArg[3] = {0.0, 0.0, 0.0};
0982   if (phiDependent)
0983     functorArg[2] = rho;
0984   else
0985     functorArg[1] = rho;
0986 
0987   for (unsigned ieta = 0; ieta < nEta; ++ieta) {
0988     const double eta(g.etaBinCenter(ieta));
0989     functorArg[0] = eta;
0990 
0991     if (phiDependent) {
0992       for (unsigned iphi = 0; iphi < nPhi; ++iphi) {
0993         functorArg[1] = g.phiBinCenter(iphi);
0994         g.uncheckedSetBin(ieta, iphi, (*f)(functorArg, 3U));
0995       }
0996     } else {
0997       const double pil = (*f)(functorArg, 2U);
0998       for (unsigned iphi = 0; iphi < nPhi; ++iphi)
0999         g.uncheckedSetBin(ieta, iphi, pil);
1000     }
1001   }
1002 }
1003 
1004 void FFTJetProducer::determinePileup() {
1005   // This function works with crisp clustering only
1006   if (!isCrisp)
1007     assert(!"Pile-up subtraction for fuzzy clustering "
1008                "is not implemented yet");
1009 
1010   // Clear the pileup vector
1011   const unsigned nJets = recoJets.size();
1012   pileup.resize(nJets);
1013   if (nJets == 0)
1014     return;
1015   const VectorLike zero;
1016   for (unsigned i = 0; i < nJets; ++i)
1017     pileup[i] = zero;
1018 
1019   // Pileup energy flow grid
1020   const fftjet::Grid2d<Real>& g(*pileupEnergyFlow);
1021   const unsigned nEta = g.nEta();
1022   const unsigned nPhi = g.nPhi();
1023   const double cellArea = g.etaBinWidth() * g.phiBinWidth();
1024 
1025   // Various calculators
1026   fftjet::Functor1<double, RecoFFTJet>& scaleCalc(*recoScaleCalcJet);
1027   fftjet::Functor1<double, RecoFFTJet>& ratioCalc(*recoScaleRatioCalcJet);
1028   fftjet::Functor1<double, RecoFFTJet>& factorCalc(*memberFactorCalcJet);
1029 
1030   // Make sure we have enough memory
1031   memFcns2dVec.resize(nJets);
1032   fftjet::AbsKernel2d** memFcns2d = &memFcns2dVec[0];
1033 
1034   doubleBuf.resize(nJets * 4U + nJets * nPhi);
1035   double* recoScales = &doubleBuf[0];
1036   double* recoScaleRatios = recoScales + nJets;
1037   double* memFactors = recoScaleRatios + nJets;
1038   double* dEta = memFactors + nJets;
1039   double* dPhiBuf = dEta + nJets;
1040 
1041   cellCountsVec.resize(nJets);
1042   unsigned* cellCounts = &cellCountsVec[0];
1043 
1044   // Go over jets and collect the necessary info
1045   for (unsigned ijet = 0; ijet < nJets; ++ijet) {
1046     const RecoFFTJet& jet(recoJets[ijet]);
1047     const fftjet::Peak& peak(jet.precluster());
1048 
1049     // Make sure we are using 2-d membership functions.
1050     // Pile-up subtraction scheme for 3-d functions should be different.
1051     fftjet::AbsMembershipFunction* m3d = dynamic_cast<fftjet::AbsMembershipFunction*>(peak.membershipFunction());
1052     if (m3d == nullptr)
1053       m3d = dynamic_cast<fftjet::AbsMembershipFunction*>(jetMembershipFunction.get());
1054     if (m3d) {
1055       assert(!"Pile-up subtraction for 3-d membership functions "
1056                    "is not implemented yet");
1057     } else {
1058       fftjet::AbsKernel2d* m2d = dynamic_cast<fftjet::AbsKernel2d*>(peak.membershipFunction());
1059       if (m2d == nullptr)
1060         m2d = dynamic_cast<fftjet::AbsKernel2d*>(jetMembershipFunction.get());
1061       assert(m2d);
1062       memFcns2d[ijet] = m2d;
1063     }
1064     recoScales[ijet] = scaleCalc(jet);
1065     recoScaleRatios[ijet] = ratioCalc(jet);
1066     memFactors[ijet] = factorCalc(jet);
1067     cellCounts[ijet] = 0U;
1068 
1069     const double jetPhi = jet.vec().Phi();
1070     for (unsigned iphi = 0; iphi < nPhi; ++iphi) {
1071       double dphi = g.phiBinCenter(iphi) - jetPhi;
1072       while (dphi > M_PI)
1073         dphi -= (2.0 * M_PI);
1074       while (dphi < -M_PI)
1075         dphi += (2.0 * M_PI);
1076       dPhiBuf[iphi * nJets + ijet] = dphi;
1077     }
1078   }
1079 
1080   // Go over all grid points and integrate
1081   // the pile-up energy density
1082   VBuilder vMaker;
1083   for (unsigned ieta = 0; ieta < nEta; ++ieta) {
1084     const double eta(g.etaBinCenter(ieta));
1085     const Real* databuf = g.data() + ieta * nPhi;
1086 
1087     // Figure out dEta for each jet
1088     for (unsigned ijet = 0; ijet < nJets; ++ijet)
1089       dEta[ijet] = eta - recoJets[ijet].vec().Eta();
1090 
1091     for (unsigned iphi = 0; iphi < nPhi; ++iphi) {
1092       double maxW(0.0);
1093       unsigned maxWJet(nJets);
1094       const double* dPhi = dPhiBuf + iphi * nJets;
1095 
1096       for (unsigned ijet = 0; ijet < nJets; ++ijet) {
1097         if (recoScaleRatios[ijet] > 0.0)
1098           memFcns2d[ijet]->setScaleRatio(recoScaleRatios[ijet]);
1099         const double f = memFactors[ijet] * (*memFcns2d[ijet])(dEta[ijet], dPhi[ijet], recoScales[ijet]);
1100         if (f > maxW) {
1101           maxW = f;
1102           maxWJet = ijet;
1103         }
1104       }
1105 
1106       if (maxWJet < nJets) {
1107         pileup[maxWJet] += vMaker(cellArea * databuf[iphi], eta, g.phiBinCenter(iphi));
1108         cellCounts[maxWJet]++;
1109       }
1110     }
1111   }
1112 }
1113 
1114 //define this as a plug-in
1115 DEFINE_FWK_MODULE(FFTJetProducer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!