Rivet/Projections/SmearedJets.hh

0001 // -*- C++ -*-
0002 #ifndef RIVET_SmearedJets_HH
0003 #define RIVET_SmearedJets_HH
0004
0005 #include "Rivet/Jet.hh"
0006 #include "Rivet/Particle.hh"
0007 #include "Rivet/Projection.hh"
0008 #include "Rivet/Projections/JetFinder.hh"
0009 #include "Rivet/Tools/SmearingFunctions.hh"
0010 #include <functional>
0011
0012 namespace Rivet {
0013
0014
0015   /// @todo Allow applying a pre-smearing cut so smearing doesn't need to be applied to below-threshold micro-jets
0016
0017
0018   /// Wrapper projection for smearing {@link Jet}s with detector resolutions and efficiencies
0019   class SmearedJets : public JetFinder {
0020   public:
0021
0022     /// @name Constructors etc.
0023     /// @{
0024
0025     /// @brief Constructor with a reco efficiency and optional tagging efficiencies
0026     ///
0027     /// @todo Add a tau-tag slot
0028     SmearedJets(const JetFinder& ja,
0029                 const JetSmearFn& smearFn,
0030                 const JetEffFn& bTagEffFn=JET_BTAG_PERFECT,
0031                 const JetEffFn& cTagEffFn=JET_CTAG_PERFECT)
0032       : SmearedJets(ja, bTagEffFn, cTagEffFn, smearFn)
0033     {    }
0034
0035
0036     /// @brief Constructor with a parameter pack of efficiency and smearing functions,
0037     /// plus optional tagging efficiencies
0038     ///
0039     /// @todo Add a tau-tag slot
0040     template <typename... Args,
0041               typename = std::enable_if_t< allArgumentsOf<JetEffSmearFn, Args...>::value >>
0042     SmearedJets(const JetFinder& ja, const JetEffFn& bTagEffFn, const JetEffFn& cTagEffFn, Args&& ... effSmearFns)
0043       : _detFns({JetEffSmearFn(std::forward<Args>(effSmearFns))...}), _bTagEffFn(bTagEffFn), _cTagEffFn(cTagEffFn)
0044     {
0045       setName("SmearedJets");
0046       declare(ja, "TruthJets");
0047       _noSmear = getEnvParam<bool>("RIVET_DISABLE_SMEARING", false);
0048     }
0049
0050     /// @todo How to include tagging effs?
0051     /// @todo Variadic eff/smear fn list?
0052     /// @todo Add a trailing Cut arg cf. SmearedParticles? -- wrap into an eff function
0053
0054
0055     /// Clone on the heap.
0056     RIVET_DEFAULT_PROJ_CLONE(SmearedJets);
0057
0058     /// @}
0059
0060     /// Import to avoid warnings about overload-hiding
0061     using Projection::operator =;
0062
0063
0064     /// Compare to another SmearedJets
0065     CmpState compare(const Projection& p) const {
0066       // Compare truth jets definitions
0067       const CmpState teq = mkPCmp(p, "TruthJets");
0068       if (teq != CmpState::EQ) return teq;
0069
0070       // Short-circuit if smearing is disabled
0071       if (_noSmear) return CmpState::EQ;
0072
0073       // Compare lists of detector functions
0074       const SmearedJets& other = dynamic_cast<const SmearedJets&>(p);
0075       const CmpState nfeq = cmp(_detFns.size(), other._detFns.size());
0076       if (nfeq != CmpState::EQ) return nfeq;
0077       for (size_t i = 0; i < _detFns.size(); ++i) {
0078         const CmpState feq = _detFns[i].cmp(other._detFns[i]);
0079         if (feq != CmpState::EQ) return feq;
0080       }
0081       return Rivet::cmp(get_address(_bTagEffFn), get_address(other._bTagEffFn)) ||
0082              Rivet::cmp(get_address(_cTagEffFn), get_address(other._cTagEffFn));
0083     }
0084
0085
0086     /// Perform the jet finding & smearing calculation
0087     void project(const Event& e) {
0088       const Jets& truthjets = apply<JetFinder>(e, "TruthJets").jetsByPt(); //truthJets();
0089
0090       // Short-circuit if smearing is disabled
0091       if (_noSmear) {
0092         _recojets = truthjets;
0093         return;
0094       }
0095
0096       // Apply jet smearing and efficiency transforms
0097       _recojets.clear(); _recojets.reserve(truthjets.size());
0098       for (const Jet& j : truthjets) {
0099         Jet jdet = j;
0100         bool keep = true;
0101         MSG_DEBUG("Truth jet: " << "mom=" << jdet.mom()/GeV << " GeV, pT=" << jdet.pT()/GeV << ", eta=" << jdet.eta());
0102         for (const JetEffSmearFn& fn : _detFns) {
0103           double jeff = -1;
0104           std::tie(jdet, jeff) = fn(jdet); // smear & eff
0105           // Re-add constituents & tags if (we assume accidentally) they were lost by the smearing function
0106           if (jdet.particles().empty() && !j.particles().empty()) jdet.particles() = j.particles();
0107           if (jdet.tags().empty() && !j.tags().empty()) jdet.tags() = j.tags();
0108           MSG_DEBUG("         ->" << "mom=" << jdet.mom()/GeV << " GeV, pT=" << jdet.pT()/GeV << ", eta=" << jdet.eta());
0109           // MSG_DEBUG("New det jet: "
0110           //           << "mom=" << jdet.mom()/GeV << " GeV, pT=" << jdet.pT()/GeV << ", eta=" << jdet.eta()
0111           //           << ", b-tag=" << boolalpha << jdet.bTagged()
0112           //           << ", c-tag=" << boolalpha << jdet.cTagged()
0113           //           << " : eff=" << 100*jeff << "%");
0114           if (jeff <= 0) { keep = false; break; } //< no need to roll expensive dice (and we deal with -ve probabilities, just in case)
0115           if (jeff < 1 && rand01() > jeff)  { keep = false; break; } //< roll dice (and deal with >1 probabilities, just in case)
0116         }
0117         if (keep) _recojets.push_back(jdet);
0118       }
0119       // Apply tagging efficiencies, using smeared kinematics as input to the tag eff functions
0120       for (Jet& j : _recojets) {
0121         // Decide whether or not there should be a b-tag on this jet
0122         const double beff = _bTagEffFn ? _bTagEffFn(j) : j.bTagged();
0123         const bool btag = beff == 1 || (beff != 0 && rand01() < beff);
0124         // Remove b-tags if needed, and add a dummy one if needed
0125         if (!btag && j.bTagged()) j.tags().erase(std::remove_if(j.tags().begin(), j.tags().end(), hasBottom), j.tags().end());
0126         if (btag && !j.bTagged()) j.tags().push_back(Particle(PID::BQUARK, j.mom())); ///< @todo Or could use the/an actual clustered b-quark momentum?
0127         // Decide whether or not there should be a c-tag on this jet
0128         const double ceff = _cTagEffFn ? _cTagEffFn(j) : j.cTagged();
0129         const bool ctag = ceff == 1 || (ceff != 0 && rand01() < beff);
0130         // Remove c-tags if needed, and add a dummy one if needed
0131         if (!ctag && j.cTagged()) j.tags().erase(std::remove_if(j.tags().begin(), j.tags().end(), hasCharm), j.tags().end());
0132         if (ctag && !j.cTagged()) j.tags().push_back(Particle(PID::CQUARK, j.mom())); ///< @todo As above... ?
0133       }
0134     }
0135
0136
0137     /// Return the full jet list for the JetFinder methods to use
0138     Jets _jets() const { return _recojets; }
0139
0140     /// Get the truth jets (sorted by pT)
0141     const Jets truthJets() const {
0142       return getProjection<JetFinder>("TruthJets").jetsByPt();
0143     }
0144
0145     /// Reset the projection. Smearing functions will be unchanged.
0146     void reset() { _recojets.clear(); }
0147
0148
0149   protected:
0150
0151     /// Smeared jets
0152     Jets _recojets;
0153
0154     /// Stored efficiency & smearing functions
0155     vector<JetEffSmearFn> _detFns;
0156
0157     /// Stored efficiency functions
0158     JetEffFn _bTagEffFn, _cTagEffFn;
0159
0160     /// Flag to disable smearing
0161     ///
0162     /// @todo Could/should be handled statically and centrally
0163     bool _noSmear{false};
0164
0165   };
0166
0167
0168 }
0169
0170 #endif