RivetInterface/src/CMS_2013_I1224539_DIJET.cc

0001 // -*- C++ -*-
0002 #include "Rivet/Analysis.hh"
0003 #include "Rivet/Projections/FinalState.hh"
0004 #include "Rivet/Projections/FastJets.hh"
0005 #include "Rivet/Projections/WFinder.hh"
0006 #include "Rivet/Projections/ZFinder.hh"
0007 #include "fastjet/tools/Filter.hh"
0008 #include "fastjet/tools/Pruner.hh"
0009
0010 namespace Rivet {
0011
0012   class CMS_2013_I1224539_DIJET : public Analysis {
0013   public:
0014     /// @name Constructors etc.
0015     //@{
0016
0017     /// Constructor
0018     CMS_2013_I1224539_DIJET()
0019         : Analysis("CMS_2013_I1224539_DIJET"),
0020           _filter(
0021               fastjet::Filter(fastjet::JetDefinition(fastjet::cambridge_algorithm, 0.3), fastjet::SelectorNHardest(3))),
0022           _trimmer(fastjet::Filter(fastjet::JetDefinition(fastjet::kt_algorithm, 0.2),
0023                                    fastjet::SelectorPtFractionMin(0.03))),
0024           _pruner(fastjet::Pruner(fastjet::cambridge_algorithm, 0.1, 0.5)) {}
0025
0026     //@}
0027
0028   public:
0029     /// @name Analysis methods
0030     //@{
0031
0032     /// Book histograms and initialise projections before the run
0033     void init() override {
0034       FinalState fs((Cuts::etaIn(-2.4, 2.4)));
0035       declare(fs, "FS");
0036
0037       // Jet collections
0038       declare(FastJets(fs, FastJets::ANTIKT, 0.7), "JetsAK7");
0039       declare(FastJets(fs, FastJets::CAM, 0.8), "JetsCA8");
0040       declare(FastJets(fs, FastJets::CAM, 1.2), "JetsCA12");
0041
0042       // Histograms
0043       for (size_t i = 0; i < N_PT_BINS_dj; ++i) {
0044         book(_h_ungroomedAvgJetMass_dj[i], i + 1 + 0 * N_PT_BINS_dj, 1, 1);
0045         book(_h_filteredAvgJetMass_dj[i], i + 1 + 1 * N_PT_BINS_dj, 1, 1);
0046         book(_h_trimmedAvgJetMass_dj[i], i + 1 + 2 * N_PT_BINS_dj, 1, 1);
0047         book(_h_prunedAvgJetMass_dj[i], i + 1 + 3 * N_PT_BINS_dj, 1, 1);
0048       }
0049     }
0050
0051     // Find the pT histogram bin index for value pt (in GeV), to hack a 2D histogram equivalent
0052     /// @todo Use a YODA axis/finder alg when available
0053     size_t findPtBin(double ptJ) {
0054       const double ptBins_dj[N_PT_BINS_dj + 1] = {220.0, 300.0, 450.0, 500.0, 600.0, 800.0, 1000.0, 1500.0};
0055       for (size_t ibin = 0; ibin < N_PT_BINS_dj; ++ibin) {
0056         if (inRange(ptJ, ptBins_dj[ibin], ptBins_dj[ibin + 1]))
0057           return ibin;
0058       }
0059       return N_PT_BINS_dj;
0060     }
0061
0062     /// Perform the per-event analysis
0063     void analyze(const Event& event) override {
0064       const double weight = 1.0;
0065
0066       // Look at events with >= 2 jets
0067       const PseudoJets& psjetsAK7 = apply<FastJets>(event, "JetsAK7").pseudoJetsByPt(50.0 * GeV);
0068       if (psjetsAK7.size() < 2)
0069         vetoEvent;
0070
0071       // Get the leading two jets and find their average pT
0072       const fastjet::PseudoJet& j0 = psjetsAK7[0];
0073       const fastjet::PseudoJet& j1 = psjetsAK7[1];
0074       double ptAvg = 0.5 * (j0.pt() + j1.pt());
0075
0076       // Find the appropriate mean pT bin and escape if needed
0077       const size_t njetBin = findPtBin(ptAvg / GeV);
0078       if (njetBin >= N_PT_BINS_dj)
0079         vetoEvent;
0080
0081       // Now run the substructure algs...
0082       fastjet::PseudoJet filtered0 = _filter(j0);
0083       fastjet::PseudoJet filtered1 = _filter(j1);
0084       fastjet::PseudoJet trimmed0 = _trimmer(j0);
0085       fastjet::PseudoJet trimmed1 = _trimmer(j1);
0086       fastjet::PseudoJet pruned0 = _pruner(j0);
0087       fastjet::PseudoJet pruned1 = _pruner(j1);
0088
0089       // ... and fill the histograms
0090       _h_ungroomedAvgJetMass_dj[njetBin]->fill(0.5 * (j0.m() + j1.m()) / GeV, weight);
0091       _h_filteredAvgJetMass_dj[njetBin]->fill(0.5 * (filtered0.m() + filtered1.m()) / GeV, weight);
0092       _h_trimmedAvgJetMass_dj[njetBin]->fill(0.5 * (trimmed0.m() + trimmed1.m()) / GeV, weight);
0093       _h_prunedAvgJetMass_dj[njetBin]->fill(0.5 * (pruned0.m() + pruned1.m()) / GeV, weight);
0094     }
0095
0096     /// Normalise histograms etc., after the run
0097     void finalize() override {
0098       const double normalizationVal = 1000;
0099       for (size_t i = 0; i < N_PT_BINS_dj; ++i) {
0100         normalize(_h_ungroomedAvgJetMass_dj[i], normalizationVal);
0101         normalize(_h_filteredAvgJetMass_dj[i], normalizationVal);
0102         normalize(_h_trimmedAvgJetMass_dj[i], normalizationVal);
0103         normalize(_h_prunedAvgJetMass_dj[i], normalizationVal);
0104       }
0105     }
0106
0107     //@}
0108
0109   private:
0110     /// @name FastJet grooming tools (configured in constructor init list)
0111     //@{
0112     const fastjet::Filter _filter;
0113     const fastjet::Filter _trimmer;
0114     const fastjet::Pruner _pruner;
0115     //@}
0116
0117     /// @name Histograms
0118     //@{
0119     enum BINS_dj {
0120       PT_220_300_dj = 0,
0121       PT_300_450_dj,
0122       PT_450_500_dj,
0123       PT_500_600_dj,
0124       PT_600_800_dj,
0125       PT_800_1000_dj,
0126       PT_1000_1500_dj,
0127       N_PT_BINS_dj
0128     };
0129     Histo1DPtr _h_ungroomedJet0pt, _h_ungroomedJet1pt;
0130     Histo1DPtr _h_ungroomedAvgJetMass_dj[N_PT_BINS_dj];
0131     Histo1DPtr _h_filteredAvgJetMass_dj[N_PT_BINS_dj];
0132     Histo1DPtr _h_trimmedAvgJetMass_dj[N_PT_BINS_dj];
0133     Histo1DPtr _h_prunedAvgJetMass_dj[N_PT_BINS_dj];
0134     //@}
0135   };
0136
0137   // The hook for the plugin system
0138   DECLARE_RIVET_PLUGIN(CMS_2013_I1224539_DIJET);
0139
0140 }  // namespace Rivet