Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:25:23 

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