Rivet/Projections/SmearedParticles.hh

0001 // -*- C++ -*-
0002 #ifndef RIVET_SmearedParticles_HH
0003 #define RIVET_SmearedParticles_HH
0004
0005 #include "Rivet/Particle.hh"
0006 #include "Rivet/Projection.hh"
0007 #include "Rivet/Projections/ParticleFinder.hh"
0008 #include "Rivet/Tools/SmearingFunctions.hh"
0009
0010 namespace Rivet {
0011
0012
0013   /// Wrapper projection for smearing {@link Jet}s with detector resolutions and efficiencies
0014   class SmearedParticles : public ParticleFinder {
0015   public:
0016
0017     /// @name Constructors etc.
0018     /// @{
0019
0020     /// @brief Constructor with a variadic ordered list of efficiency and smearing function args
0021     template<typename... Args,
0022              typename = std::enable_if_t< allArgumentsOf<ParticleEffSmearFn, Args...>::value >>
0023     SmearedParticles(const ParticleFinder& pf, Args&& ... effSmearFns)
0024       : SmearedParticles(pf, Cuts::open(), std::forward<Args>(effSmearFns) ...)
0025     {    }
0026
0027     /// @brief Constructor with a variadic ordered list of efficiency and smearing function args
0028     ///
0029     /// @note The Cut must be provided *before* the eff/smearing functions
0030     ///
0031     /// @todo Wouldn't it be nice if the Cut could also go *after* the parameter pack?
0032     template<typename... Args,
0033              typename = std::enable_if_t< allArgumentsOf<ParticleEffSmearFn, Args...>::value >>
0034     SmearedParticles(const ParticleFinder& pf, const Cut& c, Args&& ... effSmearFns)
0035       : ParticleFinder(c), _detFns({ParticleEffSmearFn(std::forward<Args>(effSmearFns))...})
0036     {
0037       setName("SmearedParticles");
0038       declare(pf, "TruthParticles");
0039       _noSmear = getEnvParam<bool>("RIVET_DISABLE_SMEARING", false);
0040     }
0041
0042
0043     /// Clone on the heap.
0044     RIVET_DEFAULT_PROJ_CLONE(SmearedParticles);
0045
0046     /// @}
0047
0048     /// Import to avoid warnings about overload-hiding
0049     using Projection::operator =;
0050
0051
0052     /// Compare to another SmearedParticles
0053     ///
0054     /// @note Comparing detector functions doesn't work for functors/lambdas,
0055     /// hence are always treated as not equivalent
0056     CmpState compare(const Projection& p) const {
0057       const SmearedParticles& other = dynamic_cast<const SmearedParticles&>(p);
0058
0059       // Compare truth particles definitions
0060       const CmpState teq = mkPCmp(other, "TruthParticles");
0061       if (teq != CmpState::EQ) return teq;
0062
0063       // Compare cuts
0064       if (_cuts != other._cuts) return CmpState::NEQ;
0065
0066       // Compare lists of detector functions
0067       if (!_noSmear) {
0068         const CmpState nfeq = cmp(_detFns.size(), other._detFns.size());
0069         MSG_TRACE("Numbers of detector functions = " << _detFns.size() << " VS " << other._detFns.size());
0070         if (nfeq != CmpState::EQ) return nfeq;
0071         for (size_t i = 0; i < _detFns.size(); ++i) {
0072           const CmpState feq = _detFns[i].cmp(other._detFns[i]);
0073           if (feq != CmpState::EQ) return feq;
0074         }
0075       }
0076
0077       // If we got this far, we're equal
0078       MSG_DEBUG("Equivalent detected! " << p.name() << ", " << this->name());
0079       return CmpState::EQ;
0080     }
0081
0082
0083     /// Perform the particle finding & smearing calculation
0084     void project(const Event& e) {
0085       const Particles& truthparticles = apply<ParticleFinder>(e, "TruthParticles").particlesByPt(); //truthParticles();
0086
0087       // Short-circuit (with cuts) if smearing is disabled
0088       if (_noSmear) {
0089         _theParticles = select(truthparticles, _cuts);
0090         return;
0091       }
0092
0093       // Apply the smearing and then reco-level cuts to each particle
0094       _theParticles.clear(); _theParticles.reserve(truthparticles.size());
0095       for (const Particle& p : truthparticles) {
0096         Particle pdet = p;
0097         double peff = -1;
0098         bool keep = true;
0099         MSG_TRACE("Number of detector functions = " << _detFns.size());
0100         for (const ParticleEffSmearFn& fn : _detFns) {
0101           std::tie(pdet, peff) = fn(pdet); // smear & eff
0102           // Test the short-circuit random numbers if possible; note handling of < 0 and > 1 probabilities
0103           if (peff <= 0 || rand01() > peff) keep = false;
0104           MSG_DEBUG("New det particle: pid=" << pdet.pid()
0105                     << ", mom=" << pdet.mom()/GeV << " GeV, "
0106                     << "pT=" << pdet.pT()/GeV << ", eta=" << pdet.eta()
0107                     << " : eff=" << 100*peff << "%, discarded=" << std::boolalpha << !keep);
0108           if (!keep) break; // discarded; no need to try more smear-eff functions
0109         }
0110         // If discarding, go straight to the next particle
0111         if (!keep) continue;
0112         // Ensure the smeared particle satisfies the cuts associated with this projection
0113         if (!_cuts->accept(pdet)) continue;
0114
0115         // Store, recording where the smearing was built from
0116         pdet.addConstituent(p); ///< @todo Is this a good idea?? What if raw particles are requested?
0117         _theParticles.push_back(pdet);
0118       }
0119     }
0120
0121     /// Get the truth particles (sorted by pT)
0122     const Particles truthParticles() const {
0123       return getProjection<ParticleFinder>("TruthParticles").particlesByPt();
0124     }
0125
0126     /// Reset the projection. Smearing functions will be unchanged.
0127     void reset() { _theParticles.clear(); }
0128
0129
0130   protected:
0131
0132     /// Stored efficiency & smearing functions
0133     vector<ParticleEffSmearFn> _detFns;
0134
0135     /// Flag to disable smearing
0136     ///
0137     /// @todo Could/should be handled statically and centrally
0138     bool _noSmear{false};
0139
0140   };
0141
0142
0143 }
0144
0145 #endif