SmearedParticles.hh

// -*- C++ -*-
#ifndef RIVET_SmearedParticles_HH
#define RIVET_SmearedParticles_HH

#include "Rivet/Particle.hh"
#include "Rivet/Projection.hh"
#include "Rivet/Projections/ParticleFinder.hh"
#include "Rivet/Tools/SmearingFunctions.hh"

namespace Rivet {


  // Recursive variadic template arg decoding
  namespace {
    template<typename T>
    vector<ParticleEffSmearFn>& toEffSmearFns(vector<ParticleEffSmearFn>& v, const T& t) {
      v.push_back(ParticleEffSmearFn(t));
      return v;
    }
    template<typename T, typename... ARGS>
    vector<ParticleEffSmearFn>& toEffSmearFns(vector<ParticleEffSmearFn>& v, const T& first, ARGS... args) {
      v.push_back(ParticleEffSmearFn(first));
      toEffSmearFns(v, args...);
      return v;
    }
  }



  class SmearedParticles : public ParticleFinder {
  public:



    SmearedParticles(const ParticleFinder& pf,
                     double eff,
                     const Cut& c=Cuts::open())
      : SmearedParticles(pf, {{eff}}, c)
    {    }

    SmearedParticles(const ParticleFinder& pf,
                     const ParticleEffFn& effFn,
                     const Cut& c=Cuts::open())
      : SmearedParticles(pf, {{effFn}}, c)
    {    }

    SmearedParticles(const ParticleFinder& pf,
                     double eff, const ParticleSmearFn& smearFn,
                     const Cut& c=Cuts::open())
      : SmearedParticles(pf, {eff, smearFn}, c)
    {    }

    SmearedParticles(const ParticleFinder& pf,
                     const ParticleSmearFn& smearFn, double eff,
                     const Cut& c=Cuts::open())
      : SmearedParticles(pf, {smearFn, eff}, c)
    {    }

    SmearedParticles(const ParticleFinder& pf,
                     const ParticleEffFn& effFn, const ParticleSmearFn& smearFn,
                     const Cut& c=Cuts::open())
      : SmearedParticles(pf, {effFn, smearFn}, c)
    {    }

    SmearedParticles(const ParticleFinder& pf,
                     const ParticleSmearFn& smearFn, const ParticleEffFn& effFn,
                     const Cut& c=Cuts::open())
      : SmearedParticles(pf, {smearFn, effFn}, c)
    {    }

    SmearedParticles(const ParticleFinder& pf,
                     const vector<ParticleEffSmearFn>& effSmearFns,
                     const Cut& c=Cuts::open())
      : ParticleFinder(c),
        _detFns(effSmearFns)
    {
      setName("SmearedParticles");
      addProjection(pf, "TruthParticles");
    }

    SmearedParticles(const ParticleFinder& pf,
                     const initializer_list<ParticleEffSmearFn>& effSmearFns,
                     const Cut& c=Cuts::open())
      : SmearedParticles(pf, vector<ParticleEffSmearFn>{effSmearFns}, c)
    {    }

    template<typename... ARGS>
    SmearedParticles(const ParticleFinder& pf, const Cut& c, ARGS... effSmearFns)
      : SmearedParticles(pf, toEffSmearFns(_detFns, effSmearFns...), c)
    {    }


    DEFAULT_RIVET_PROJ_CLONE(SmearedParticles);



    int compare(const Projection& p) const {
      const SmearedParticles& other = dynamic_cast<const SmearedParticles&>(p);

      // Compare truth particles definitions
      const int teq = mkPCmp(other, "TruthParticles");
      if (teq != EQUIVALENT) return teq;

      // Compare lists of detector functions
      const int nfeq = cmp(_detFns.size(), other._detFns.size());
      if (nfeq != EQUIVALENT) return nfeq;
      for (size_t i = 0; i < _detFns.size(); ++i) {
        const int feq = _detFns[i].cmp(other._detFns[i]);
        if (feq != EQUIVALENT) return feq;
      }

      // If we got this far, we're equal
      return EQUIVALENT;
    }


    void project(const Event& e) {
      // Copying and filtering
      const Particles& truthparticles = apply<ParticleFinder>(e, "TruthParticles").particlesByPt();
      _theParticles.clear(); _theParticles.reserve(truthparticles.size());
      for (const Particle& p : truthparticles) {
        Particle pdet = p;
        double peff = -1;
        bool keep = true;
        for (const ParticleEffSmearFn& fn : _detFns) {
          tie(pdet, peff) = fn(pdet); // smear & eff
          MSG_DEBUG("New det particle: pid=" << pdet.pid()
                    << ", mom=" << pdet.mom()/GeV << " GeV, "
                    << "pT=" << pdet.pT()/GeV << ", eta=" << pdet.eta()
                    << " : eff=" << 100*peff << "%");
          if (peff <= 0) { keep = false; break; } //< no need to roll expensive dice (and we deal with -ve probabilities, just in case)
          if (peff < 1 && rand01() > peff)  { keep = false; break; } //< roll dice (and deal with >1 probabilities, just in case)
        }
        if (keep) {
          pdet.addConstituent(p); //< record where the smearing was built from
          _theParticles.push_back(pdet);
        }
      }
    }


    void reset() { _theParticles.clear(); }


  private:

    vector<ParticleEffSmearFn> _detFns;

  };


}

#endif