FFTJetProducers/plugins/FFTJetDijetFilter.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    FFTJetProducers
0004 // Class:      FFTJetDijetFilter
0005 //
0006 /**\class FFTJetDijetFilter FFTJetDijetFilter.cc RecoJets/FFTJetProducers/plugins/FFTJetDijetFilter.cc
0007
0008  Description: selects good dijet events looking only at the clustering tree
0009
0010  Implementation:
0011      [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  Igor Volobouev
0015 //         Created:  Thu Jul 18 19:19:40 CDT 2012
0016 //
0017 //
0018 #include <cmath>
0019 #include <cassert>
0020 #include <iostream>
0021
0022 #include "fftjet/EquidistantSequence.hh"
0023 #include "fftjet/ProximityClusteringTree.hh"
0024 #include "fftjet/SparseClusteringTree.hh"
0025 #include "fftjet/PeakEtaPhiDistance.hh"
0026 #include "fftjet/peakEtLifetime.hh"
0027
0028 // framework include files
0029 #include "FWCore/Framework/interface/Frameworkfwd.h"
0030 #include "FWCore/Framework/interface/stream/EDFilter.h"
0031 #include "FWCore/Framework/interface/Event.h"
0032 #include "FWCore/Framework/interface/MakerMacros.h"
0033 #include "FWCore/Utilities/interface/Exception.h"
0034 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0035
0036 #include "RecoJets/FFTJetAlgorithms/interface/fftjetTypedefs.h"
0037 #include "RecoJets/FFTJetAlgorithms/interface/clusteringTreeConverters.h"
0038 #include "RecoJets/FFTJetProducers/interface/FFTJetParameterParser.h"
0039
0040 #include "DataFormats/Math/interface/deltaPhi.h"
0041 #include "DataFormats/Common/interface/View.h"
0042 #include "DataFormats/Common/interface/Handle.h"
0043
0044 #define init_param(type, varname) varname(ps.getParameter<type>(#varname))
0045
0046 using namespace fftjetcms;
0047
0048 //
0049 // class declaration
0050 //
0051 class FFTJetDijetFilter : public edm::stream::EDFilter<> {
0052 public:
0053   typedef fftjet::ProximityClusteringTree<fftjet::Peak, long> ClusteringTree;
0054   typedef fftjet::SparseClusteringTree<fftjet::Peak, long> SparseTree;
0055
0056   explicit FFTJetDijetFilter(const edm::ParameterSet&);
0057   FFTJetDijetFilter() = delete;
0058   FFTJetDijetFilter(const FFTJetDijetFilter&) = delete;
0059   FFTJetDijetFilter& operator=(const FFTJetDijetFilter&) = delete;
0060   ~FFTJetDijetFilter() override;
0061
0062 private:
0063   typedef reco::PattRecoTree<float, reco::PattRecoPeak<float> > StoredTree;
0064
0065   bool filter(edm::Event& iEvent, const edm::EventSetup& iSetup) override;
0066
0067   template <class Ptr>
0068   inline void checkConfig(const Ptr& ptr, const char* message) const {
0069     if (ptr.get() == nullptr)
0070       throw cms::Exception("FFTJetBadConfig") << message << std::endl;
0071   }
0072
0073   inline double peakPt(const fftjet::Peak& peak) const {
0074     const double s = peak.scale();
0075     return ptConversionFactor * s * s * peak.magnitude();
0076   }
0077
0078   // Module parameters
0079   edm::InputTag treeLabel;
0080   double ptConversionFactor;
0081   double fixedScale;
0082   double completeEventScale;
0083   double min1to0PtRatio;
0084   double minDeltaPhi;
0085   double maxThirdJetFraction;
0086   double minPt0;
0087   double minPt1;
0088   double maxPeakEta;
0089   bool insertCompleteEvent;
0090
0091   // Distance calculator for the clustering tree
0092   std::unique_ptr<fftjet::AbsDistanceCalculator<fftjet::Peak> > distanceCalc;
0093
0094   // Scales used
0095   std::unique_ptr<std::vector<double> > iniScales;
0096
0097   // The clustering trees
0098   ClusteringTree* clusteringTree;
0099   SparseTree* sparseTree;
0100
0101   // Space for peaks
0102   std::vector<fftjet::Peak> peaks;
0103
0104   // Space for sparse tree nodes
0105   std::vector<unsigned> nodes;
0106
0107   // pass/fail decision counters
0108   unsigned long nPassed;
0109   unsigned long nRejected;
0110
0111   edm::EDGetTokenT<StoredTree> treeToken;
0112 };
0113
0114 //
0115 // constructors and destructor
0116 //
0117 FFTJetDijetFilter::FFTJetDijetFilter(const edm::ParameterSet& ps)
0118     : init_param(edm::InputTag, treeLabel),
0119       init_param(double, ptConversionFactor),
0120       init_param(double, fixedScale),
0121       init_param(double, completeEventScale),
0122       init_param(double, min1to0PtRatio),
0123       init_param(double, minDeltaPhi),
0124       init_param(double, maxThirdJetFraction),
0125       init_param(double, minPt0),
0126       init_param(double, minPt1),
0127       init_param(double, maxPeakEta),
0128       init_param(bool, insertCompleteEvent),
0129       clusteringTree(nullptr),
0130       sparseTree(nullptr),
0131       nPassed(0),
0132       nRejected(0) {
0133   // Parse the set of scales
0134   iniScales = fftjet_ScaleSet_parser(ps.getParameter<edm::ParameterSet>("InitialScales"));
0135   checkConfig(iniScales, "invalid set of scales");
0136   std::sort(iniScales->begin(), iniScales->end(), std::greater<double>());
0137
0138   // Parse the distance calculator
0139   const edm::ParameterSet& TreeDistanceCalculator(ps.getParameter<edm::ParameterSet>("TreeDistanceCalculator"));
0140   distanceCalc = fftjet_DistanceCalculator_parser(TreeDistanceCalculator);
0141   checkConfig(distanceCalc, "invalid tree distance calculator");
0142
0143   // Create the clustering tree
0144   clusteringTree = new ClusteringTree(distanceCalc.get());
0145   sparseTree = new SparseTree();
0146
0147   treeToken = consumes<StoredTree>(treeLabel);
0148 }
0149
0150 FFTJetDijetFilter::~FFTJetDijetFilter() {
0151   delete sparseTree;
0152   delete clusteringTree;
0153 }
0154
0155 // ------------ method called to produce the data  ------------
0156 bool FFTJetDijetFilter::filter(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0157   edm::Handle<StoredTree> input;
0158   iEvent.getByToken(treeToken, input);
0159
0160   // Convert the stored tree into a normal FFTJet clustering tree
0161   // and extract the set of peaks at the requested scale
0162   const double eventScale = insertCompleteEvent ? completeEventScale : 0.0;
0163   if (input->isSparse()) {
0164     sparsePeakTreeFromStorable(*input, iniScales.get(), eventScale, sparseTree);
0165     sparseTree->sortNodes();
0166     fftjet::updateSplitMergeTimes(*sparseTree, sparseTree->minScale(), sparseTree->maxScale());
0167     const unsigned usedLevel = sparseTree->getLevel(fixedScale);
0168     sparseTree->getLevelNodes(usedLevel, &nodes);
0169     const unsigned numNodes = nodes.size();
0170     peaks.clear();
0171     peaks.reserve(numNodes);
0172     for (unsigned i = 0; i < numNodes; ++i)
0173       peaks.push_back(sparseTree->uncheckedNode(nodes[i]).getCluster());
0174   } else {
0175     densePeakTreeFromStorable(*input, iniScales.get(), eventScale, clusteringTree);
0176     const unsigned usedLevel = clusteringTree->getLevel(fixedScale);
0177     double actualScale = 0.0;
0178     long dummyInfo;
0179     clusteringTree->getLevelData(usedLevel, &actualScale, &peaks, &dummyInfo);
0180   }
0181
0182   // Get out if we don't have two clusters
0183   const unsigned nClusters = peaks.size();
0184   if (nClusters < 2) {
0185     ++nRejected;
0186     return false;
0187   }
0188   std::sort(peaks.begin(), peaks.end(), std::greater<fftjet::Peak>());
0189
0190   // Calculate all quantities needed to make the pass/fail decision
0191   const double pt0 = peakPt(peaks[0]);
0192   const double pt1 = peakPt(peaks[1]);
0193   const double dphi = reco::deltaPhi(peaks[0].phi(), peaks[1].phi());
0194   const double ptratio = pt1 / pt0;
0195   double thirdJetFraction = 0.0;
0196   if (nClusters > 2)
0197     thirdJetFraction = peakPt(peaks[2]) / (pt0 + pt1);
0198
0199   // Calculate the cut
0200   const bool pass = pt0 > minPt0 && pt1 > minPt1 && ptratio > min1to0PtRatio && std::abs(dphi) > minDeltaPhi &&
0201                     thirdJetFraction < maxThirdJetFraction && std::abs(peaks[0].eta()) < maxPeakEta &&
0202                     std::abs(peaks[1].eta()) < maxPeakEta;
0203   if (pass)
0204     ++nPassed;
0205   else
0206     ++nRejected;
0207   return pass;
0208 }
0209
0210 //define this as a plug-in
0211 DEFINE_FWK_MODULE(FFTJetDijetFilter);