File indexing completed on 2022-05-21 03:43:06
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
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
0028 #include "RecoJets/FFTJetProducers/plugins/FFTJetProducer.h"
0029
0030
0031 #include "FWCore/Framework/interface/MakerMacros.h"
0032 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0033
0034
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
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
0052
0053
0054
0055
0056
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 }
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
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
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
0170
0171
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
0185
0186
0187
0188 }
0189
0190 FFTJetProducer::~FFTJetProducer() {}
0191
0192
0193
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
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& ,
0213 edm::Event& iEvent,
0214 const edm::EventSetup& ,
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
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
0251
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
0265
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
0303
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
0312
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
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
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
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
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
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
0452
0453
0454
0455 converged = checkConvergence(recoJets, iterJets);
0456
0457
0458 iterJets.swap(recoJets);
0459 nJets = recoJets.size();
0460 }
0461
0462
0463 if (preclusters.size() != nJets) {
0464 assert(nJets < preclusters.size());
0465 removeFakePreclusters();
0466 assert(preclusters.size() == nJets);
0467 }
0468
0469
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
0477
0478
0479
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
0515
0516
0517
0518 const unsigned mask = clusterMask[i];
0519 assert(mask <= nJets);
0520 constituents[mask].push_back(inputCollection->ptrAt(candIdx[i]));
0521 }
0522 }
0523
0524
0525
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
0534
0535
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
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
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
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
0575
0576
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
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
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
0604 const double pt = jet.pt();
0605 if (pt > previousPt)
0606 sorted = false;
0607 previousPt = pt;
0608 }
0609
0610
0611 if (!sorted)
0612 std::sort(jets->begin(), jets->end(), LocalSortByPt());
0613
0614
0615 iEvent.put(std::move(jets), outputLabel);
0616 }
0617
0618 void FFTJetProducer::saveResults(edm::Event& ev, const edm::EventSetup& iSetup, const unsigned nPreclustersFound) {
0619
0620 jet_type_switch(writeJets, ev, iSetup);
0621
0622
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
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
0654 void FFTJetProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0655
0656 if (storeInSinglePrecision())
0657 loadSparseTreeData<float>(iEvent, input_recotree_token_f_);
0658 else
0659 loadSparseTreeData<double>(iEvent, input_recotree_token_d_);
0660
0661
0662 if (assignConstituents || !(useGriddedAlgorithm && reuseExistingGrid))
0663 loadInputCollection(iEvent);
0664
0665
0666 if (useGriddedAlgorithm) {
0667 if (reuseExistingGrid) {
0668 if (loadEnergyFlow(iEvent, energyFlow))
0669 buildGridAlg();
0670 } else
0671 discretizeEnergyFlow();
0672 }
0673
0674
0675 sparseTree.occupancyInScaleSpace(*peakSelector, &occupancy);
0676
0677
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
0689 prepareRecombinationScales();
0690
0691
0692
0693
0694 assignMembershipFunctions(&preclusters);
0695
0696
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
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
0713 if (maxIterations > 1U && !recoJets.empty()) {
0714
0715
0716
0717
0718
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
0729
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
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
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
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
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
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
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
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
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
0816 void FFTJetProducer::beginStream(edm::StreamID) {
0817 const edm::ParameterSet& ps(myConfiguration);
0818
0819
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
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
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
0862 recoScaleCalcPeak = parse_recoScaleCalcPeak(ps);
0863 checkConfig(recoScaleCalcPeak,
0864 "invalid spec for the "
0865 "reconstruction scale calculator from peaks");
0866
0867
0868 recoScaleRatioCalcPeak = parse_recoScaleRatioCalcPeak(ps);
0869 checkConfig(recoScaleRatioCalcPeak,
0870 "invalid spec for the "
0871 "reconstruction scale ratio calculator from peaks");
0872
0873
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
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
0905
0906
0907
0908
0909
0910
0911 std::vector<int> matchTable;
0912 const unsigned nmatched = matchOneToOne(recoJets, preclusters, JetToPeakDistance(), &matchTable);
0913
0914
0915
0916 assert(nmatched == recoJets.size());
0917
0918
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
1006 if (!isCrisp)
1007 assert(!"Pile-up subtraction for fuzzy clustering "
1008 "is not implemented yet");
1009
1010
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
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
1026 fftjet::Functor1<double, RecoFFTJet>& scaleCalc(*recoScaleCalcJet);
1027 fftjet::Functor1<double, RecoFFTJet>& ratioCalc(*recoScaleRatioCalcJet);
1028 fftjet::Functor1<double, RecoFFTJet>& factorCalc(*memberFactorCalcJet);
1029
1030
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
1045 for (unsigned ijet = 0; ijet < nJets; ++ijet) {
1046 const RecoFFTJet& jet(recoJets[ijet]);
1047 const fftjet::Peak& peak(jet.precluster());
1048
1049
1050
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
1081
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
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
1115 DEFINE_FWK_MODULE(FFTJetProducer);