File indexing completed on 2024-04-06 12:25:23
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 esLoader_.acquireToken(pileupTableRecord, consumesCollector());
0185
0186
0187
0188
0189
0190 }
0191
0192 FFTJetProducer::~FFTJetProducer() {}
0193
0194
0195
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
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& ,
0215 edm::Event& iEvent,
0216 const edm::EventSetup& ,
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
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
0253
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
0267
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
0305
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
0314
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
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
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
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
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
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
0454
0455
0456
0457 converged = checkConvergence(recoJets, iterJets);
0458
0459
0460 iterJets.swap(recoJets);
0461 nJets = recoJets.size();
0462 }
0463
0464
0465 if (preclusters.size() != nJets) {
0466 assert(nJets < preclusters.size());
0467 removeFakePreclusters();
0468 assert(preclusters.size() == nJets);
0469 }
0470
0471
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
0479
0480
0481
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
0517
0518
0519
0520 const unsigned mask = clusterMask[i];
0521 assert(mask <= nJets);
0522 constituents[mask].push_back(inputCollection->ptrAt(candIdx[i]));
0523 }
0524 }
0525
0526
0527
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
0536
0537
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
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
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
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
0577
0578
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
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
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
0606 const double pt = jet.pt();
0607 if (pt > previousPt)
0608 sorted = false;
0609 previousPt = pt;
0610 }
0611
0612
0613 if (!sorted)
0614 std::sort(jets->begin(), jets->end(), LocalSortByPt());
0615
0616
0617 iEvent.put(std::move(jets), outputLabel);
0618 }
0619
0620 void FFTJetProducer::saveResults(edm::Event& ev, const edm::EventSetup& iSetup, const unsigned nPreclustersFound) {
0621
0622 jet_type_switch(writeJets, ev, iSetup);
0623
0624
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
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
0656 void FFTJetProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0657
0658 if (storeInSinglePrecision())
0659 loadSparseTreeData<float>(iEvent, input_recotree_token_f_);
0660 else
0661 loadSparseTreeData<double>(iEvent, input_recotree_token_d_);
0662
0663
0664 if (assignConstituents || !(useGriddedAlgorithm && reuseExistingGrid))
0665 loadInputCollection(iEvent);
0666
0667
0668 if (useGriddedAlgorithm) {
0669 if (reuseExistingGrid) {
0670 if (loadEnergyFlow(iEvent, energyFlow))
0671 buildGridAlg();
0672 } else
0673 discretizeEnergyFlow();
0674 }
0675
0676
0677 sparseTree.occupancyInScaleSpace(*peakSelector, &occupancy);
0678
0679
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
0691 prepareRecombinationScales();
0692
0693
0694
0695
0696 assignMembershipFunctions(&preclusters);
0697
0698
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
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
0715 if (maxIterations > 1U && !recoJets.empty()) {
0716
0717
0718
0719
0720
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
0731
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
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
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
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
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
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
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
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
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
0818 void FFTJetProducer::beginStream(edm::StreamID) {
0819 const edm::ParameterSet& ps(myConfiguration);
0820
0821
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
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
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
0864 recoScaleCalcPeak = parse_recoScaleCalcPeak(ps);
0865 checkConfig(recoScaleCalcPeak,
0866 "invalid spec for the "
0867 "reconstruction scale calculator from peaks");
0868
0869
0870 recoScaleRatioCalcPeak = parse_recoScaleRatioCalcPeak(ps);
0871 checkConfig(recoScaleRatioCalcPeak,
0872 "invalid spec for the "
0873 "reconstruction scale ratio calculator from peaks");
0874
0875
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
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
0907
0908
0909
0910
0911
0912
0913 std::vector<int> matchTable;
0914 const unsigned nmatched = matchOneToOne(recoJets, preclusters, JetToPeakDistance(), &matchTable);
0915
0916
0917
0918 assert(nmatched == recoJets.size());
0919
0920
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
1007 if (!isCrisp)
1008 assert(!"Pile-up subtraction for fuzzy clustering "
1009 "is not implemented yet");
1010
1011
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
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
1027 fftjet::Functor1<double, RecoFFTJet>& scaleCalc(*recoScaleCalcJet);
1028 fftjet::Functor1<double, RecoFFTJet>& ratioCalc(*recoScaleRatioCalcJet);
1029 fftjet::Functor1<double, RecoFFTJet>& factorCalc(*memberFactorCalcJet);
1030
1031
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
1046 for (unsigned ijet = 0; ijet < nJets; ++ijet) {
1047 const RecoFFTJet& jet(recoJets[ijet]);
1048 const fftjet::Peak& peak(jet.precluster());
1049
1050
1051
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
1082
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
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
1116 DEFINE_FWK_MODULE(FFTJetProducer);