File indexing completed on 2023-03-17 11:18:18
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019 #include <fstream>
0020 #include <memory>
0021
0022
0023 #include "fftjet/ProximityClusteringTree.hh"
0024 #include "fftjet/ClusteringSequencer.hh"
0025 #include "fftjet/ClusteringTreeSparsifier.hh"
0026 #include "fftjet/FrequencyKernelConvolver.hh"
0027 #include "fftjet/FrequencySequentialConvolver.hh"
0028 #include "fftjet/DiscreteGauss1d.hh"
0029 #include "fftjet/DiscreteGauss2d.hh"
0030
0031
0032 #include "FWCore/Framework/interface/Frameworkfwd.h"
0033 #include "FWCore/Framework/interface/MakerMacros.h"
0034 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0035
0036 #include "DataFormats/Common/interface/View.h"
0037 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
0038 #include "DataFormats/JetReco/interface/GenJetCollection.h"
0039 #include "DataFormats/JetReco/interface/PFJetCollection.h"
0040 #include "DataFormats/JetReco/interface/BasicJetCollection.h"
0041 #include "DataFormats/JetReco/interface/TrackJetCollection.h"
0042 #include "DataFormats/Candidate/interface/Candidate.h"
0043
0044
0045 #include "DataFormats/JetReco/interface/DiscretizedEnergyFlow.h"
0046
0047
0048 #include "RecoJets/FFTJetProducers/interface/FFTJetParameterParser.h"
0049
0050
0051 #include "RecoJets/FFTJetAlgorithms/interface/clusteringTreeConverters.h"
0052
0053
0054 #include "RecoJets/FFTJetAlgorithms/interface/gridConverters.h"
0055
0056
0057 #include "RecoJets/FFTJetProducers/interface/FFTJetInterface.h"
0058
0059 using namespace fftjetcms;
0060
0061
0062
0063
0064 class FFTJetPatRecoProducer : public FFTJetInterface {
0065 public:
0066 explicit FFTJetPatRecoProducer(const edm::ParameterSet&);
0067 FFTJetPatRecoProducer() = delete;
0068 FFTJetPatRecoProducer(const FFTJetPatRecoProducer&) = delete;
0069 FFTJetPatRecoProducer& operator=(const FFTJetPatRecoProducer&) = delete;
0070 ~FFTJetPatRecoProducer() override;
0071
0072 protected:
0073
0074 typedef fftjet::ProximityClusteringTree<fftjet::Peak, long> ClusteringTree;
0075 typedef fftjet::SparseClusteringTree<fftjet::Peak, long> SparseTree;
0076 typedef fftjet::ClusteringSequencer<Real> Sequencer;
0077 typedef fftjet::ClusteringTreeSparsifier<fftjet::Peak, long> Sparsifier;
0078
0079
0080 void produce(edm::Event&, const edm::EventSetup&) override;
0081
0082 void buildKernelConvolver(const edm::ParameterSet&);
0083 fftjet::PeakFinder buildPeakFinder(const edm::ParameterSet&);
0084
0085 template <class Real>
0086 void buildSparseProduct(edm::Event&) const;
0087
0088 template <class Real>
0089 void buildDenseProduct(edm::Event&) const;
0090
0091
0092 ClusteringTree* clusteringTree;
0093
0094
0095
0096
0097 std::unique_ptr<Sequencer> sequencer;
0098 std::unique_ptr<Sparsifier> sparsifier;
0099
0100
0101 std::unique_ptr<MyFFTEngine> engine;
0102 std::unique_ptr<MyFFTEngine> anotherEngine;
0103
0104
0105 std::unique_ptr<fftjet::AbsFrequencyKernel> kernel2d;
0106 std::unique_ptr<fftjet::AbsFrequencyKernel1d> etaKernel;
0107 std::unique_ptr<fftjet::AbsFrequencyKernel1d> phiKernel;
0108
0109
0110 std::unique_ptr<fftjet::AbsConvolverBase<Real> > convolver;
0111
0112
0113 std::unique_ptr<fftjet::Functor1<bool, fftjet::Peak> > peakSelector;
0114
0115
0116 std::unique_ptr<fftjet::AbsDistanceCalculator<fftjet::Peak> > distanceCalc;
0117
0118
0119 SparseTree sparseTree;
0120
0121
0122
0123
0124 const double completeEventDataCutoff;
0125
0126
0127 const bool makeClusteringTree;
0128
0129
0130
0131 const bool verifyDataConversion;
0132
0133
0134 const bool storeDiscretizationGrid;
0135
0136
0137 const bool sparsify;
0138
0139 private:
0140
0141 std::ofstream externalGridStream;
0142 bool storeGridsExternally;
0143 fftjet::Grid2d<float>* extGrid;
0144 };
0145
0146
0147
0148
0149 FFTJetPatRecoProducer::FFTJetPatRecoProducer(const edm::ParameterSet& ps)
0150 : FFTJetInterface(ps),
0151 clusteringTree(nullptr),
0152 completeEventDataCutoff(ps.getParameter<double>("completeEventDataCutoff")),
0153 makeClusteringTree(ps.getParameter<bool>("makeClusteringTree")),
0154 verifyDataConversion(ps.getUntrackedParameter<bool>("verifyDataConversion", false)),
0155 storeDiscretizationGrid(ps.getParameter<bool>("storeDiscretizationGrid")),
0156 sparsify(ps.getParameter<bool>("sparsify")),
0157 extGrid(nullptr) {
0158
0159 if (makeClusteringTree) {
0160 if (storeInSinglePrecision())
0161 produces<reco::PattRecoTree<float, reco::PattRecoPeak<float> > >(outputLabel);
0162 else
0163 produces<reco::PattRecoTree<double, reco::PattRecoPeak<double> > >(outputLabel);
0164 }
0165 if (storeDiscretizationGrid)
0166 produces<reco::DiscretizedEnergyFlow>(outputLabel);
0167
0168
0169 const std::string externalGridFile(ps.getParameter<std::string>("externalGridFile"));
0170 storeGridsExternally = !externalGridFile.empty();
0171 if (storeGridsExternally) {
0172 externalGridStream.open(externalGridFile.c_str(), std::ios_base::out | std::ios_base::binary);
0173 if (!externalGridStream.is_open())
0174 throw cms::Exception("FFTJetBadConfig")
0175 << "FFTJetPatRecoProducer failed to open file " << externalGridFile << std::endl;
0176 }
0177
0178 if (!makeClusteringTree && !storeDiscretizationGrid && !storeGridsExternally) {
0179 throw cms::Exception("FFTJetBadConfig")
0180 << "FFTJetPatRecoProducer is not configured to produce anything" << std::endl;
0181 }
0182
0183
0184
0185
0186 energyFlow = fftjet_Grid2d_parser(ps.getParameter<edm::ParameterSet>("GridConfiguration"));
0187 checkConfig(energyFlow, "invalid discretization grid");
0188
0189
0190
0191 buildKernelConvolver(ps);
0192
0193
0194 peakSelector = fftjet_PeakSelector_parser(ps.getParameter<edm::ParameterSet>("PeakSelectorConfiguration"));
0195 checkConfig(peakSelector, "invalid peak selector");
0196
0197
0198 std::unique_ptr<std::vector<double> > iniScales =
0199 fftjet_ScaleSet_parser(ps.getParameter<edm::ParameterSet>("InitialScales"));
0200 checkConfig(iniScales, "invalid set of scales");
0201
0202
0203 const unsigned maxAdaptiveScales = ps.getParameter<unsigned>("maxAdaptiveScales");
0204 const double minAdaptiveRatioLog = ps.getParameter<double>("minAdaptiveRatioLog");
0205 if (minAdaptiveRatioLog <= 0.0)
0206 throw cms::Exception("FFTJetBadConfig") << "bad adaptive ratio logarithm limit" << std::endl;
0207
0208
0209
0210 if (getEventScale() > 0.0) {
0211 const double cs = getEventScale();
0212 const unsigned nscales = iniScales->size();
0213 for (unsigned i = 0; i < nscales; ++i)
0214 if (cs >= (*iniScales)[i])
0215 throw cms::Exception("FFTJetBadConfig") << "incompatible scale for complete event" << std::endl;
0216 }
0217
0218
0219 sequencer = std::make_unique<Sequencer>(
0220 convolver.get(), peakSelector.get(), buildPeakFinder(ps), *iniScales, maxAdaptiveScales, minAdaptiveRatioLog);
0221
0222
0223 const edm::ParameterSet& SparsifierConfiguration(ps.getParameter<edm::ParameterSet>("SparsifierConfiguration"));
0224 sparsifier = fftjet_ClusteringTreeSparsifier_parser(SparsifierConfiguration);
0225 checkConfig(sparsifier, "invalid sparsifier parameters");
0226
0227
0228 const edm::ParameterSet& TreeDistanceCalculator(ps.getParameter<edm::ParameterSet>("TreeDistanceCalculator"));
0229 distanceCalc = fftjet_DistanceCalculator_parser(TreeDistanceCalculator);
0230 checkConfig(distanceCalc, "invalid tree distance calculator");
0231
0232
0233 clusteringTree = new ClusteringTree(distanceCalc.get());
0234 }
0235
0236 void FFTJetPatRecoProducer::buildKernelConvolver(const edm::ParameterSet& ps) {
0237
0238
0239 const std::vector<double> etaDependentScaleFactors(ps.getParameter<std::vector<double> >("etaDependentScaleFactors"));
0240
0241
0242 const bool use2dKernel = etaDependentScaleFactors.empty();
0243 if (!use2dKernel)
0244 if (etaDependentScaleFactors.size() != energyFlow->nEta())
0245 throw cms::Exception("FFTJetBadConfig") << "invalid number of eta-dependent scale factors" << std::endl;
0246
0247
0248 double kernelEtaScale = ps.getParameter<double>("kernelEtaScale");
0249 const double kernelPhiScale = ps.getParameter<double>("kernelPhiScale");
0250 if (kernelEtaScale <= 0.0 || kernelPhiScale <= 0.0)
0251 throw cms::Exception("FFTJetBadConfig") << "invalid kernel scale" << std::endl;
0252
0253
0254
0255 kernelEtaScale *= (2.0 * M_PI / (energyFlow->etaMax() - energyFlow->etaMin()));
0256
0257
0258 const bool fixEfficiency = ps.getParameter<bool>("fixEfficiency");
0259
0260
0261 unsigned convolverMinBin = 0, convolverMaxBin = 0;
0262 if (fixEfficiency || !use2dKernel) {
0263 convolverMinBin = ps.getParameter<unsigned>("convolverMinBin");
0264 convolverMaxBin = ps.getParameter<unsigned>("convolverMaxBin");
0265 }
0266
0267 if (use2dKernel) {
0268
0269 engine = std::make_unique<MyFFTEngine>(energyFlow->nEta(), energyFlow->nPhi());
0270
0271
0272 kernel2d = std::unique_ptr<fftjet::AbsFrequencyKernel>(
0273 new fftjet::DiscreteGauss2d(kernelEtaScale, kernelPhiScale, energyFlow->nEta(), energyFlow->nPhi()));
0274
0275
0276 convolver = std::unique_ptr<fftjet::AbsConvolverBase<Real> >(new fftjet::FrequencyKernelConvolver<Real, Complex>(
0277 engine.get(), kernel2d.get(), convolverMinBin, convolverMaxBin));
0278 } else {
0279
0280 engine = std::make_unique<MyFFTEngine>(1, energyFlow->nEta());
0281 anotherEngine = std::make_unique<MyFFTEngine>(1, energyFlow->nPhi());
0282
0283
0284 etaKernel =
0285 std::unique_ptr<fftjet::AbsFrequencyKernel1d>(new fftjet::DiscreteGauss1d(kernelEtaScale, energyFlow->nEta()));
0286
0287 phiKernel =
0288 std::unique_ptr<fftjet::AbsFrequencyKernel1d>(new fftjet::DiscreteGauss1d(kernelPhiScale, energyFlow->nPhi()));
0289
0290
0291 convolver = std::unique_ptr<fftjet::AbsConvolverBase<Real> >(
0292 new fftjet::FrequencySequentialConvolver<Real, Complex>(engine.get(),
0293 anotherEngine.get(),
0294 etaKernel.get(),
0295 phiKernel.get(),
0296 etaDependentScaleFactors,
0297 convolverMinBin,
0298 convolverMaxBin,
0299 fixEfficiency));
0300 }
0301 }
0302
0303 fftjet::PeakFinder FFTJetPatRecoProducer::buildPeakFinder(const edm::ParameterSet& ps) {
0304 const double peakFinderMaxEta = ps.getParameter<double>("peakFinderMaxEta");
0305 if (peakFinderMaxEta <= 0.0)
0306 throw cms::Exception("FFTJetBadConfig") << "invalid peak finder eta cut" << std::endl;
0307 const double maxMagnitude = ps.getParameter<double>("peakFinderMaxMagnitude");
0308 int minBin = energyFlow->getEtaBin(-peakFinderMaxEta);
0309 if (minBin < 0)
0310 minBin = 0;
0311 int maxBin = energyFlow->getEtaBin(peakFinderMaxEta) + 1;
0312 if (maxBin < 0)
0313 maxBin = 0;
0314 return fftjet::PeakFinder(maxMagnitude, true, minBin, maxBin);
0315 }
0316
0317 FFTJetPatRecoProducer::~FFTJetPatRecoProducer() {
0318
0319
0320 if (storeGridsExternally)
0321 externalGridStream.close();
0322 delete clusteringTree;
0323 delete extGrid;
0324 }
0325
0326
0327
0328
0329 template <class Real>
0330 void FFTJetPatRecoProducer::buildSparseProduct(edm::Event& ev) const {
0331 typedef reco::PattRecoTree<Real, reco::PattRecoPeak<Real> > StoredTree;
0332
0333 auto tree = std::make_unique<StoredTree>();
0334
0335 sparsePeakTreeToStorable(sparseTree, sequencer->maxAdaptiveScales(), tree.get());
0336
0337
0338 if (verifyDataConversion && !storeInSinglePrecision()) {
0339 SparseTree check;
0340 const std::vector<double>& scalesUsed(sequencer->getInitialScales());
0341 sparsePeakTreeFromStorable(*tree, &scalesUsed, getEventScale(), &check);
0342 if (sparseTree != check)
0343 throw cms::Exception("FFTJetInterface") << "Data conversion failed for sparse clustering tree" << std::endl;
0344 }
0345
0346 ev.put(std::move(tree), outputLabel);
0347 }
0348
0349 template <class Real>
0350 void FFTJetPatRecoProducer::buildDenseProduct(edm::Event& ev) const {
0351 typedef reco::PattRecoTree<Real, reco::PattRecoPeak<Real> > StoredTree;
0352
0353 auto tree = std::make_unique<StoredTree>();
0354
0355 densePeakTreeToStorable(*clusteringTree, sequencer->maxAdaptiveScales(), tree.get());
0356
0357
0358 if (verifyDataConversion && !storeInSinglePrecision()) {
0359 ClusteringTree check(distanceCalc.get());
0360 const std::vector<double>& scalesUsed(sequencer->getInitialScales());
0361 densePeakTreeFromStorable(*tree, &scalesUsed, getEventScale(), &check);
0362 if (*clusteringTree != check)
0363 throw cms::Exception("FFTJetInterface") << "Data conversion failed for dense clustering tree" << std::endl;
0364 }
0365
0366 ev.put(std::move(tree), outputLabel);
0367 }
0368
0369
0370 void FFTJetPatRecoProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0371 loadInputCollection(iEvent);
0372 discretizeEnergyFlow();
0373
0374 if (makeClusteringTree) {
0375 sequencer->run(*energyFlow, clusteringTree);
0376 if (getEventScale() > 0.0)
0377 sequencer->insertCompleteEvent(getEventScale(), *energyFlow, clusteringTree, completeEventDataCutoff);
0378
0379 if (sparsify) {
0380 sparsifier->sparsify(*clusteringTree, &sparseTree);
0381
0382
0383
0384
0385
0386
0387 if (storeInSinglePrecision())
0388 buildSparseProduct<float>(iEvent);
0389 else
0390 buildSparseProduct<double>(iEvent);
0391 } else {
0392 if (storeInSinglePrecision())
0393 buildDenseProduct<float>(iEvent);
0394 else
0395 buildDenseProduct<double>(iEvent);
0396 }
0397 }
0398
0399 if (storeDiscretizationGrid) {
0400 const fftjet::Grid2d<Real>& g(*energyFlow);
0401
0402 auto flow = std::make_unique<reco::DiscretizedEnergyFlow>(
0403 g.data(), g.title(), g.etaMin(), g.etaMax(), g.phiBin0Edge(), g.nEta(), g.nPhi());
0404
0405 if (verifyDataConversion) {
0406 fftjet::Grid2d<Real> check(
0407 flow->nEtaBins(), flow->etaMin(), flow->etaMax(), flow->nPhiBins(), flow->phiBin0Edge(), flow->title());
0408 check.blockSet(flow->data(), flow->nEtaBins(), flow->nPhiBins());
0409 assert(g == check);
0410 }
0411
0412 iEvent.put(std::move(flow), outputLabel);
0413 }
0414
0415 if (storeGridsExternally) {
0416 if (extGrid)
0417 copy_Grid2d_data(extGrid, *energyFlow);
0418 else
0419 extGrid = convert_Grid2d_to_float(*energyFlow);
0420 if (!extGrid->write(externalGridStream)) {
0421 throw cms::Exception("FFTJetPatRecoProducer::produce")
0422 << "Failed to write grid data into an external file" << std::endl;
0423 }
0424 }
0425 }
0426
0427
0428 DEFINE_FWK_MODULE(FFTJetPatRecoProducer);