PercentileProjection.hh

// -*- C++ -*-
#ifndef RIVET_PERCENTILEPROJECTION_HH
#define RIVET_PERCENTILEPROJECTION_HH
 
#include "Rivet/Projections/SingleValueProjection.hh"
#include "Rivet/Tools/RivetYODA.hh"
#include <map>
 
namespace Rivet {
 
  enum class PercentileOrder { INCREASING, DECREASING };
 

  class PercentileProjection : public SingleValueProjection {
  public:
 
    using SingleValueProjection::operator=;

    PercentileProjection(const SingleValueProjection& sv,
                         const Histo1D& calhist,
                         PercentileOrder pctorder = PercentileOrder::DECREASING)
        : _calhist("EMPTY"), _increasing(pctorder == PercentileOrder::INCREASING) {
      setName("PercentileProjection");
      declare(sv, "OBSERVABLE");
      //if ( !calhist ) return;
      MSG_DEBUG("Constructing PercentileProjection from " << calhist.path());
      _calhist = calhist.path();
      int N = calhist.numBins();
      double sum = calhist.sumW();
      if (_increasing) {
        double acc = 0.0;
        for (int i = 0; i <= N; ++i) {
          acc += calhist.bin(i).sumW();
          _table.insert(make_pair(calhist.bin(i).xMax(), 100.0 * acc / sum));
        }
      }
      else {
        double acc = 0.0;
        for (int i = N + 1; i > 0; --i) {
          acc += calhist.bin(i).sumW();
          _table.insert(make_pair(calhist.bin(i).xMin(), 100.0 * acc / sum));
        }
      }
      if (getLog().isActive(Log::DEBUG)) {
        MSG_DEBUG("Mapping from observable to percentile:");
        for (auto p : _table) {
          std::cout << std::setw(16) << p.first << " -> " << std::setw(16) << p.second << "%" << std::endl;
          if (not _increasing and p.second <= 0) break;
          if (_increasing and p.second >= 100) break;
        }
      }
    }
 
 
    // Constructor taking a SingleValueProjection and a calibration
    // histogram. If increasing it means that low values corresponds to
    // lower percentiles.
    PercentileProjection(const SingleValueProjection& sv,
                         const Estimate1D& calest,
                         PercentileOrder pctorder = PercentileOrder::DECREASING)
        : _calhist("EMPTY"), _increasing(pctorder == PercentileOrder::INCREASING) {
      declare(sv, "OBSERVABLE");
 
      //if ( !calest ) return;
      MSG_DEBUG("Constructing PercentileProjection from " << calest.path());
      _calhist = calest.path();
      int N = calest.numBins();
      double sum = 0.0;
      for (const auto& b : calest.bins()) sum += b.val();
 
      double acc = 0.0;
      if (_increasing) {
        _table.insert(make_pair(calest.bin(1).xMin(), 100.0 * acc / sum));
        for (int i = 0; i < N; ++i) {
          acc += calest.bin(i + 1).val();
          _table.insert(make_pair(calest.bin(i + 1).xMax(), 100.0 * acc / sum));
        }
      }
      else {
        _table.insert(make_pair(calest.bin(N).xMax(), 100.0 * acc / sum));
        for (int i = N - 1; i >= 0; --i) {
          acc += calest.bin(i + 1).val();
          _table.insert(make_pair(calest.bin(i + 1).xMin(), 100.0 * acc / sum));
        }
      }
    }
 
 
    RIVET_DEFAULT_PROJ_CLONE(PercentileProjection);

    using Projection::operator=;
 
 
    // The projection function takes the assigned SingeValueProjection
    // and sets the value of this projection to the corresponding
    // percentile. If no calibration has been provided, -1 will be
    // returned. If values are outside of the calibration histogram, 0
    // or 100 will be returned.
    void project(const Event& e) {
      clear();
      if (_table.empty()) return;
      auto& pobs = apply<SingleValueProjection>(e, "OBSERVABLE");
      double obs = pobs();
      double pcnt = lookup(obs);
      if (pcnt >= 0.0) setValue(pcnt);
      MSG_DEBUG("Observable(" << pobs.name() << ")=" << std::setw(16) << obs
                              << "-> Percentile=" << std::setw(16) << pcnt << "%");
    }
 
    // Standard comparison function.
    CmpState compare(const Projection& p) const {
      const PercentileProjection pp = dynamic_cast<const PercentileProjection&>(p);
      return mkNamedPCmp(p, "OBSERVABLE") || cmp(_increasing, pp._increasing) || cmp(_calhist, pp._calhist);
    }
 
 
  protected:
 
    // The (interpolated) lookup table
    double lookup(double obs) const {
      auto low = _table.upper_bound(obs);
      if (low == _table.end()) return _increasing ? 100.0 : 0.0;
      if (low == _table.begin()) return _increasing ? 0.0 : 100.0;
      auto high = low--;
      return low->second + (obs - low->first) * (high->second - low->second) / (high->first - low->first);
    }
 
    // Astring identifying the calibration histogram.
    string _calhist;
 
    // A lookup table to find (by interpolation) the percentile given
    // the value of the underlying SingleValueProjection.
    map<double, double> _table;
 
    // A flag to say whether the distribution should be integrated from
    // below or above.
    bool _increasing;
  };
 
 
}
 
#endif