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     /// @note The Cut must be provided *before* the eff/smearing functions
0029     /// @todo Wouldn't it be nice if the Cut could also go *after* the parameter pack?
0030     template<typename... Args,
0031              typename = std::enable_if_t< allArgumentsOf<ParticleEffSmearFn, Args...>::value >>
0032     SmearedParticles(const ParticleFinder& pf, const Cut& c, Args&& ... effSmearFns)
0033       : ParticleFinder(c), _detFns({ParticleEffSmearFn(std::forward<Args>(effSmearFns))...})
0034     {
0035       setName("SmearedParticles");
0036       declare(pf, "TruthParticles");
0037     }
0038
0039
0040     /// Clone on the heap.
0041     RIVET_DEFAULT_PROJ_CLONE(SmearedParticles);
0042
0043     /// @}
0044
0045     /// Import to avoid warnings about overload-hiding
0046     using Projection::operator =;
0047
0048
0049     /// Compare to another SmearedParticles
0050     ///
0051     /// @note Comparing detector functions doesn't work for functors/lambdas,
0052     /// hence are always treated as not equivalent
0053     CmpState compare(const Projection& p) const {
0054       const SmearedParticles& other = dynamic_cast<const SmearedParticles&>(p);
0055
0056       // Compare truth particles definitions
0057       const CmpState teq = mkPCmp(other, "TruthParticles");
0058       if (teq != CmpState::EQ) return teq;
0059
0060       // Compare cuts
0061       if (_cuts != other._cuts) return CmpState::NEQ;
0062
0063       // Compare lists of detector functions
0064       const CmpState nfeq = cmp(_detFns.size(), other._detFns.size());
0065       MSG_TRACE("Numbers of detector functions = " << _detFns.size() << " VS " << other._detFns.size());
0066       if (nfeq != CmpState::EQ) return nfeq;
0067       for (size_t i = 0; i < _detFns.size(); ++i) {
0068         const CmpState feq = _detFns[i].cmp(other._detFns[i]);
0069         if (feq != CmpState::EQ) return feq;
0070       }
0071
0072       // If we got this far, we're equal
0073       MSG_DEBUG("Equivalent detected! " << p.name() << ", " << this->name());
0074       return CmpState::EQ;
0075     }
0076
0077
0078     /// Perform the particle finding & smearing calculation
0079     void project(const Event& e) {
0080       // Copying and filtering
0081       const Particles& truthparticles = apply<ParticleFinder>(e, "TruthParticles").particlesByPt(); //truthParticles();
0082       _theParticles.clear(); _theParticles.reserve(truthparticles.size());
0083       for (const Particle& p : truthparticles) {
0084         Particle pdet = p;
0085         double peff = -1;
0086         bool keep = true;
0087         MSG_TRACE("Number of detector functions = " << _detFns.size());
0088         for (const ParticleEffSmearFn& fn : _detFns) {
0089           std::tie(pdet, peff) = fn(pdet); // smear & eff
0090           // Test the short-circuit random numbers if possible; note handling of < 0 and > 1 probabilities
0091           if (peff <= 0 || rand01() > peff) keep = false;
0092           MSG_DEBUG("New det particle: pid=" << pdet.pid()
0093                     << ", mom=" << pdet.mom()/GeV << " GeV, "
0094                     << "pT=" << pdet.pT()/GeV << ", eta=" << pdet.eta()
0095                     << " : eff=" << 100*peff << "%, discarded=" << std::boolalpha << !keep);
0096           if (!keep) break; // discarded; no need to try more smear-eff functions
0097         }
0098         // If discarding, go straight to the next particle
0099         if (!keep) continue;
0100         //Ensure the smeared particle satisfies the cuts associated with this projection
0101         if (!_cuts->accept(pdet)) continue;
0102
0103         // Store, recording where the smearing was built from
0104         pdet.addConstituent(p); ///< @todo Is this a good idea?? What if raw particles are requested?
0105         _theParticles.push_back(pdet);
0106       }
0107     }
0108
0109     /// Get the truth particles (sorted by pT)
0110     const Particles truthParticles() const {
0111       return getProjection<ParticleFinder>("TruthParticles").particlesByPt();
0112     }
0113
0114     /// Reset the projection. Smearing functions will be unchanged.
0115     void reset() { _theParticles.clear(); }
0116
0117
0118   protected:
0119
0120     /// Stored efficiency & smearing functions
0121     vector<ParticleEffSmearFn> _detFns;
0122
0123   };
0124
0125
0126 }
0127
0128 #endif