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     }
0048
0049     /// @todo How to include tagging effs?
0050     /// @todo Variadic eff/smear fn list?
0051     /// @todo Add a trailing Cut arg cf. SmearedParticles? -- wrap into an eff function
0052
0053
0054     /// Clone on the heap.
0055     RIVET_DEFAULT_PROJ_CLONE(SmearedJets);
0056
0057     /// @}
0058
0059     /// Import to avoid warnings about overload-hiding
0060     using Projection::operator =;
0061
0062
0063     /// Compare to another SmearedJets
0064     CmpState compare(const Projection& p) const {
0065       // Compare truth jets definitions
0066       const CmpState teq = mkPCmp(p, "TruthJets");
0067       if (teq != CmpState::EQ) return teq;
0068
0069       // Compare lists of detector functions
0070       const SmearedJets& other = dynamic_cast<const SmearedJets&>(p);
0071       const CmpState nfeq = cmp(_detFns.size(), other._detFns.size());
0072       if (nfeq != CmpState::EQ) return nfeq;
0073       for (size_t i = 0; i < _detFns.size(); ++i) {
0074         const CmpState feq = _detFns[i].cmp(other._detFns[i]);
0075         if (feq != CmpState::EQ) return feq;
0076       }
0077       return Rivet::cmp(get_address(_bTagEffFn), get_address(other._bTagEffFn)) ||
0078              Rivet::cmp(get_address(_cTagEffFn), get_address(other._cTagEffFn));
0079     }
0080
0081
0082     /// Perform the jet finding & smearing calculation
0083     void project(const Event& e) {
0084       // Copying and filtering
0085       const Jets& truthjets = apply<JetFinder>(e, "TruthJets").jetsByPt(); //truthJets();
0086       _recojets.clear(); _recojets.reserve(truthjets.size());
0087       // Apply jet smearing and efficiency transforms
0088       for (const Jet& j : truthjets) {
0089         Jet jdet = j;
0090         bool keep = true;
0091         MSG_DEBUG("Truth jet: " << "mom=" << jdet.mom()/GeV << " GeV, pT=" << jdet.pT()/GeV << ", eta=" << jdet.eta());
0092         for (const JetEffSmearFn& fn : _detFns) {
0093           double jeff = -1;
0094           std::tie(jdet, jeff) = fn(jdet); // smear & eff
0095           // Re-add constituents & tags if (we assume accidentally) they were lost by the smearing function
0096           if (jdet.particles().empty() && !j.particles().empty()) jdet.particles() = j.particles();
0097           if (jdet.tags().empty() && !j.tags().empty()) jdet.tags() = j.tags();
0098           MSG_DEBUG("         ->" << "mom=" << jdet.mom()/GeV << " GeV, pT=" << jdet.pT()/GeV << ", eta=" << jdet.eta());
0099           // MSG_DEBUG("New det jet: "
0100           //           << "mom=" << jdet.mom()/GeV << " GeV, pT=" << jdet.pT()/GeV << ", eta=" << jdet.eta()
0101           //           << ", b-tag=" << boolalpha << jdet.bTagged()
0102           //           << ", c-tag=" << boolalpha << jdet.cTagged()
0103           //           << " : eff=" << 100*jeff << "%");
0104           if (jeff <= 0) { keep = false; break; } //< no need to roll expensive dice (and we deal with -ve probabilities, just in case)
0105           if (jeff < 1 && rand01() > jeff)  { keep = false; break; } //< roll dice (and deal with >1 probabilities, just in case)
0106         }
0107         if (keep) _recojets.push_back(jdet);
0108       }
0109       // Apply tagging efficiencies, using smeared kinematics as input to the tag eff functions
0110       for (Jet& j : _recojets) {
0111         // Decide whether or not there should be a b-tag on this jet
0112         const double beff = _bTagEffFn ? _bTagEffFn(j) : j.bTagged();
0113         const bool btag = beff == 1 || (beff != 0 && rand01() < beff);
0114         // Remove b-tags if needed, and add a dummy one if needed
0115         if (!btag && j.bTagged()) j.tags().erase(std::remove_if(j.tags().begin(), j.tags().end(), hasBottom), j.tags().end());
0116         if (btag && !j.bTagged()) j.tags().push_back(Particle(PID::BQUARK, j.mom())); ///< @todo Or could use the/an actual clustered b-quark momentum?
0117         // Decide whether or not there should be a c-tag on this jet
0118         const double ceff = _cTagEffFn ? _cTagEffFn(j) : j.cTagged();
0119         const bool ctag = ceff == 1 || (ceff != 0 && rand01() < beff);
0120         // Remove c-tags if needed, and add a dummy one if needed
0121         if (!ctag && j.cTagged()) j.tags().erase(std::remove_if(j.tags().begin(), j.tags().end(), hasCharm), j.tags().end());
0122         if (ctag && !j.cTagged()) j.tags().push_back(Particle(PID::CQUARK, j.mom())); ///< @todo As above... ?
0123       }
0124     }
0125
0126
0127     /// Return the full jet list for the JetFinder methods to use
0128     Jets _jets() const { return _recojets; }
0129
0130     /// Get the truth jets (sorted by pT)
0131     const Jets truthJets() const {
0132       return getProjection<JetFinder>("TruthJets").jetsByPt();
0133     }
0134
0135     /// Reset the projection. Smearing functions will be unchanged.
0136     void reset() { _recojets.clear(); }
0137
0138
0139   protected:
0140
0141     /// Smeared jets
0142     Jets _recojets;
0143
0144     /// Stored efficiency & smearing functions
0145     vector<JetEffSmearFn> _detFns;
0146
0147     /// Stored efficiency functions
0148     JetEffFn _bTagEffFn, _cTagEffFn;
0149
0150   };
0151
0152
0153 }
0154
0155 #endif