IsoDeposit.h

Go to the documentation of this file.

#ifndef RecoCandidate_IsoDeposit_H
#define RecoCandidate_IsoDeposit_H

#include "DataFormats/RecoCandidate/interface/IsoDepositDirection.h"
#include "DataFormats/Math/interface/Vector3D.h"
#include "FWCore/Utilities/interface/Exception.h"
#include "FWCore/Utilities/interface/thread_safety_macros.h"
#include <map>
#include <cmath>
#include <string>
#include <vector>
#include <typeinfo>
#include <atomic>

namespace reco {
  namespace isodeposit {
    struct AbsVeto {
      virtual ~AbsVeto() {}
      virtual bool veto(double eta, double phi, float value) const = 0;
      virtual void centerOn(double eta, double phi) {
        throw cms::Exception("Not Implemented") << "This AbsVeto implementation (" << typeid(this).name()
                                                << ") does not support the centerOn(eta,phi) method";
      }
    };
    typedef std::vector<AbsVeto*> AbsVetos;
  }  // namespace isodeposit
}  // namespace reco

namespace reco {

  class IsoDeposit {
  public:
    typedef isodeposit::Direction Direction;
    typedef isodeposit::AbsVeto AbsVeto;
    typedef isodeposit::AbsVetos AbsVetos;
    typedef Direction::Distance Distance;
    typedef std::multimap<Distance, float> DepositsMultimap;
    typedef DepositsMultimap::const_iterator DepIterator;

    // old style vetos
    struct Veto {
      Direction vetoDir;
      float dR;
      Veto() {}
      Veto(Direction dir, double d) : vetoDir(dir), dR(d) {}
    };
    typedef std::vector<Veto> Vetos;

    IsoDeposit(double eta = 0, double phi = 0);
    IsoDeposit(const Direction& candDirection);

    virtual ~IsoDeposit(){};

    const Direction& direction() const { return theDirection; }
    double eta() const { return theDirection.eta(); }
    double phi() const { return theDirection.phi(); }

    const Veto& veto() const { return theVeto; }
    void setVeto(const Veto& aVeto) { theVeto = aVeto; }

    void addDeposit(double dr, double deposit);  // FIXME - temporary for backward compatibility
    void addDeposit(const Direction& depDir, double deposit);

    double depositWithin(double coneSize,               //dR in which deposit is computed
                         const Vetos& vetos = Vetos(),  //additional vetos
                         bool skipDepositVeto = false   //skip exclusion of veto
                         ) const;

    double depositWithin(Direction dir,
                         double coneSize,               //dR in which deposit is computed
                         const Vetos& vetos = Vetos(),  //additional vetos
                         bool skipDepositVeto = false   //skip exclusion of veto
                         ) const;

    std::pair<double, int> depositAndCountWithin(double coneSize,               //dR in which deposit is computed
                                                 const Vetos& vetos = Vetos(),  //additional vetos
                                                 double threshold = -1e+36,     //threshold on counted deposits
                                                 bool skipDepositVeto = false   //skip exclusion of veto
                                                 ) const;

    std::pair<double, int> depositAndCountWithin(Direction dir,                 //wrt another direction
                                                 double coneSize,               //dR in which deposit is computed
                                                 const Vetos& vetos = Vetos(),  //additional vetos
                                                 double threshold = -1e+36,     //threshold on deposits
                                                 bool skipDepositVeto = false   //skip exclusion of veto
                                                 ) const;

    double depositWithin(double coneSize,              //dR in which deposit is computed
                         const AbsVetos& vetos,        //additional vetos
                         bool skipDepositVeto = false  //skip exclusion of veto
                         ) const;

    std::pair<double, int> depositAndCountWithin(double coneSize,              //dR in which deposit is computed
                                                 const AbsVetos& vetos,        //additional vetos
                                                 bool skipDepositVeto = false  //skip exclusion of veto
                                                 ) const;

    double candEnergy() const { return theCandTag; }

    void addCandEnergy(double et) { theCandTag += et; }

    std::string print() const;

    class const_iterator {
    public:
      const const_iterator& operator++() {
        ++it_;
        cacheReady_ = false;
        return *this;
      }
      const const_iterator* operator->() const { return this; }
      float dR() const { return it_->first.deltaR; }
      float eta() const {
        if (!cacheReady_)
          doDir();
        return cache_.eta();
      }
      float phi() const {
        if (!cacheReady_)
          doDir();
        return cache_.phi();
      }
      float value() const { return it_->second; }
      bool operator!=(const const_iterator& it2) { return it2.it_ != it_; }
      friend class IsoDeposit;

    private:
      typedef Direction::Distance Distance;
      void doDir() const {
        if (!cacheReady_)
          cache_ = parent_->direction() + it_->first;
        cacheReady_ = true;
      }
      const_iterator(const IsoDeposit* parent, std::multimap<Distance, float>::const_iterator it)
          : parent_(parent), it_(it), cache_(), cacheReady_(false) {}
      const reco::IsoDeposit* parent_;
      std::multimap<Distance, float>::const_iterator it_;
      CMS_THREAD_SAFE mutable Direction cache_;
      mutable std::atomic<bool> cacheReady_;
    };
    const_iterator begin() const { return const_iterator(this, theDeposits.begin()); }
    const_iterator end() const { return const_iterator(this, theDeposits.end()); }

    class SumAlgo {
    public:
      SumAlgo() : sum_(0) {}
      void operator+=(DepIterator deposit) { sum_ += deposit->second; }
      void operator+=(double deposit) { sum_ += deposit; }
      double result() const { return sum_; }

    private:
      double sum_;
    };
    class CountAlgo {
    public:
      CountAlgo() : count_(0) {}
      void operator+=(DepIterator deposit) { count_++; }
      void operator+=(double deposit) { count_++; }
      double result() const { return count_; }

    private:
      size_t count_;
    };
    class Sum2Algo {
    public:
      Sum2Algo() : sum2_(0) {}
      void operator+=(DepIterator deposit) { sum2_ += deposit->second * deposit->second; }
      void operator+=(double deposit) { sum2_ += deposit * deposit; }
      double result() const { return sum2_; }

    private:
      double sum2_;
    };
    class MaxAlgo {
    public:
      MaxAlgo() : max_(0) {}
      void operator+=(DepIterator deposit) {
        if (deposit->second > max_)
          max_ = deposit->second;
      }
      void operator+=(double deposit) {
        if (deposit > max_)
          max_ = deposit;
      }
      double result() const { return max_; }

    private:
      double max_;
    };

    class MeanDRAlgo {
    public:
      MeanDRAlgo() : sum_(0.), count_(0.) {}
      void operator+=(DepIterator deposit) {
        sum_ += deposit->first.deltaR;
        count_ += 1.0;
      }
      double result() const { return sum_ / std::max(1., count_); }

    private:
      double sum_;
      double count_;
    };

    class SumDRAlgo {
    public:
      SumDRAlgo() : sum_(0.) {}
      void operator+=(DepIterator deposit) { sum_ += deposit->first.deltaR; }
      double result() const { return sum_; }

    private:
      double sum_;
    };

    template <typename Algo>
    double algoWithin(double coneSize,                     //dR in which deposit is computed
                      const AbsVetos& vetos = AbsVetos(),  //additional vetos
                      bool skipDepositVeto = false         //skip exclusion of veto
                      ) const;
    template <typename Algo>
    double algoWithin(const Direction&,
                      double coneSize,                     //dR in which deposit is computed
                      const AbsVetos& vetos = AbsVetos(),  //additional vetos
                      bool skipDepositVeto = false         //skip exclusion of veto
                      ) const;
    // count of the non-vetoed deposits in the cone
    double countWithin(double coneSize,                     //dR in which deposit is computed
                       const AbsVetos& vetos = AbsVetos(),  //additional vetos
                       bool skipDepositVeto = false         //skip exclusion of veto
                       ) const;
    // sum of the non-vetoed deposits in the cone
    double sumWithin(double coneSize,                     //dR in which deposit is computed
                     const AbsVetos& vetos = AbsVetos(),  //additional vetos
                     bool skipDepositVeto = false         //skip exclusion of veto
                     ) const;
    // sum of the non-vetoed deposits in the cone w.r.t. other direction
    double sumWithin(const Direction& dir,
                     double coneSize,                      //dR in which deposit is computed
                     const AbsVetos& vetos = AbsVetos(),   //additional vetos
                     bool skipDepositVeto = false          //skip exclusion of veto
                     ) const;                              // sum of the squares of the non-vetoed deposits in the cone
    double sum2Within(double coneSize,                     //dR in which deposit is computed
                      const AbsVetos& vetos = AbsVetos(),  //additional vetos
                      bool skipDepositVeto = false         //skip exclusion of veto
                      ) const;
    // maximum value among the non-vetoed deposits in the cone
    double maxWithin(double coneSize,                     //dR in which deposit is computed
                     const AbsVetos& vetos = AbsVetos(),  //additional vetos
                     bool skipDepositVeto = false         //skip exclusion of veto
                     ) const;
    // maximum value among the non-vetoed deposits in the cone
    double nearestDR(double coneSize,                     //dR in which deposit is computed
                     const AbsVetos& vetos = AbsVetos(),  //additional vetos
                     bool skipDepositVeto = false         //skip exclusion of veto
                     ) const;

  private:
    Direction theDirection;

    Veto theVeto;

    float theCandTag;


    DepositsMultimap theDeposits;
  };

}  // namespace reco

template <typename Algo>
double reco::IsoDeposit::algoWithin(double coneSize, const AbsVetos& vetos, bool skipDepositVeto) const {
  using namespace reco::isodeposit;
  Algo algo;
  typedef AbsVetos::const_iterator IV;
  IV ivEnd = vetos.end();

  Distance maxDistance = {float(coneSize), 999.f};
  typedef DepositsMultimap::const_iterator IM;
  IM imLoc = theDeposits.upper_bound(maxDistance);
  for (IM im = theDeposits.begin(); im != imLoc; ++im) {
    bool vetoed = false;
    Direction dirDep = theDirection + im->first;
    for (IV iv = vetos.begin(); iv < ivEnd; ++iv) {
      if ((*iv)->veto(dirDep.eta(), dirDep.phi(), im->second)) {
        vetoed = true;
        break;
      }
    }
    if (!vetoed) {
      if (skipDepositVeto || (dirDep.deltaR(theVeto.vetoDir) > theVeto.dR)) {
        algo += im;
      }
    }
  }
  return algo.result();
}

template <typename Algo>
double reco::IsoDeposit::algoWithin(const Direction& dir,
                                    double coneSize,
                                    const AbsVetos& vetos,
                                    bool skipDepositVeto) const {
  using namespace reco::isodeposit;
  Algo algo;
  typedef AbsVetos::const_iterator IV;
  IV ivEnd = vetos.end();
  typedef DepositsMultimap::const_iterator IM;
  IM imLoc = theDeposits.end();
  for (IM im = theDeposits.begin(); im != imLoc; ++im) {
    bool vetoed = false;
    Direction dirDep = theDirection + im->first;
    Distance newDist = dirDep - dir;
    if (newDist.deltaR > coneSize)
      continue;
    for (IV iv = vetos.begin(); iv < ivEnd; ++iv) {
      if ((*iv)->veto(dirDep.eta(), dirDep.phi(), im->second)) {
        vetoed = true;
        break;
      }
    }
    if (!vetoed) {
      if (skipDepositVeto || (dirDep.deltaR(theVeto.vetoDir) > theVeto.dR)) {
        algo += im;
      }
    }
  }
  return algo.result();
}

#endif