ConstIntCut.h

Go to the documentation of this file.

// \author Darko Veberic
// \date 2014

#ifndef _cic_ConstIntCut_h_
#define _cic_ConstIntCut_h_

#include <cic/Event.h>
#include <utl/Accumulator.h>
#include <utl/Bin.h>
#include <vector>
#include <algorithm>
#include <functional>


namespace cic {

  struct ConstIntCut {

    typedef std::function<double(double)> Attenuation;
    typedef std::function<double(double, double)> SFunction;
    typedef utl::Accumulator::Collect<Event> BinType;

    ConstIntCut(const std::vector<Event>& e,
                const SFunction& sFunc,
                const double maxSin2Theta,
                const unsigned int nBinsSin2Theta,
                const double sRange) :
      fMaxSin2Theta(maxSin2Theta),
      fSin2ThetaBins(
        utl::EquidistantInNorm<Event>(
          0,
          maxSin2Theta,
          nBinsSin2Theta,
          [](const Event& e) { return e.fSin2Theta; }
        )
      ),
      fSFunc(sFunc),
      fNEvents(e.size())
    {
      fSin2ThetaBins.Apply(e);
      for (BinType& b : fSin2ThetaBins.GetAllBins())
        b.Sort(
          [&sRange](const Event& a, const Event& b) {
            return a.fLnS1000 + sRange*a.fLnS1000Err > b.fLnS1000 + sRange*b.fLnS1000Err;
          }
        );
    }

    unsigned int GetNBins() const { return fSin2ThetaBins.GetNBins(); }

  private:
    static
    void
    SortByFirst(std::vector<std::pair<double, double>>& pairs)
    {
      typedef std::pair<double, double> DD;
      std::sort(pairs.begin(), pairs.end(),
                [](const DD& a, const DD& b) { return a.first < b.first; });
    }

    void
    ResetCounters()
    {
      fNEventsAboveCut = 0;
      fSin2ThetaMinMaxN =
        utl::Accumulator::MinMax<double>(std::numeric_limits<double>::infinity(),
                                         -std::numeric_limits<double>::infinity());
      fLnSMinMaxN = fSin2ThetaMinMaxN;
    }

    double
    GetMinLnAttenuation(const Attenuation& att,
                        const double mins2, const double maxs2)
    {
      using namespace std;

      const unsigned int ns2 = 10;
      const double ds2 = (maxs2 - mins2) / ns2;
      utl::Accumulator::Min<double> minLnAtt(log(att(mins2)));
      for (unsigned int i = 1; i <= ns2; ++i)
        minLnAtt(log(att(mins2 + i*ds2)));
      const double eps = 0.01;
      return minLnAtt.GetMin() - eps;
    }

  public:
    std::vector<double>
    MakePDF(const Attenuation& att, const double lnSCut)
    {
      using namespace std;

      const unsigned int n = fSin2ThetaBins.GetNBins();
      vector<double> pdf;
      pdf.reserve(n);
      ResetCounters();
      for (unsigned int i = 0; i < n; ++i) {
        const double minLnAtt =
          GetMinLnAttenuation(att, fSin2ThetaBins.GetBinLowerEdge(i),
                                   fSin2ThetaBins.GetBinUpperEdge(i));
        const BinType& b = fSin2ThetaBins[i];
        double sum = 0;
        for (const Event& e : b.GetCollection()) {
          const double cut1000 = lnSCut + log(att(e.fSin2Theta));
          const double w = fSFunc(e.fLnS1000 - cut1000, e.fLnS1000Err);
          if (!w && !fSFunc(e.fLnS1000 - (lnSCut + minLnAtt), e.fLnS1000Err))
            break;
          sum += w;
        }
        pdf.push_back(sum);
        fNEventsAboveCut += sum;
        fSin2ThetaMinMaxN(sum);
      }
      return pdf;
    }

    template<class Binning>
    std::vector<std::vector<double>>
    MakePDF2D(const Attenuation& att, const Binning& lnS38Binning)
    {
      using namespace std;
      using namespace utl;
      using namespace utl::Accumulator;

      Bin<Sum, Binning> binLnS38(lnS38Binning);
      const unsigned int nlns = binLnS38.GetNBins();
      const double minLnSCut = binLnS38.GetBinLowerEdge(0);
      const unsigned int ns2 = fSin2ThetaBins.GetNBins();
      vector<vector<double>> ret(nlns, vector<double>(ns2));
      const vector<Collect<Event>>& s2Bins = fSin2ThetaBins.GetAllBins();
      for (unsigned int is2 = 0; is2 < ns2; ++is2) {
        const double minLnAtt =
          GetMinLnAttenuation(att, fSin2ThetaBins.GetBinLowerEdge(is2),
                                   fSin2ThetaBins.GetBinUpperEdge(is2));
        for (const Event& e : s2Bins[is2].GetCollection()) {
          const double lnS38 = e.fLnS1000 - log(att(e.fSin2Theta));
          const double err = e.fLnS1000Err;
          double wup = fSFunc(lnS38 - binLnS38.GetBinUpperEdge(nlns-1), err);
          for (int ilns = nlns-1; ilns >= 0; --ilns) {
            const double wdn =
              fSFunc(lnS38 - binLnS38.GetBinLowerEdge(ilns), err);
            const double w = wdn - wup;
            if (w)
              binLnS38[ilns](w);
            wup = wdn;
          }
          if (!wup && !fSFunc(e.fLnS1000 - (minLnSCut + minLnAtt), err))
            break;
        }
        for (unsigned int ilns = 0; ilns < nlns; ++ilns)
          ret[ilns][is2] = binLnS38[ilns].GetSum();
        binLnS38.Clear();
      }

      ResetCounters();
      for (const auto& row : ret) {
        double rowSum = 0;
        for (const double n : row) {
          rowSum += n;
          fSin2ThetaMinMaxN(n);
        }
        fNEventsAboveCut += rowSum;
        fLnSMinMaxN(rowSum);
      }

      return ret;
    }

    std::vector<std::pair<double, double>>
    MakeCDF(const Attenuation& att, const double lnSCut)
    {
      using namespace std;

      const unsigned int n = fSin2ThetaBins.GetNBins();
      vector<pair<double, double>> wcdf;
      ResetCounters();
      for (unsigned int i = 0; i < n; ++i) {
        const double minLnAtt =
          GetMinLnAttenuation(att, fSin2ThetaBins.GetBinLowerEdge(i),
                                   fSin2ThetaBins.GetBinUpperEdge(i));
        for (const Event& e : fSin2ThetaBins[i].GetCollection()) {
          const double s2 = e.fSin2Theta;
          const double cut1000 = lnSCut + log(att(s2));
          const double w = fSFunc(e.fLnS1000 - cut1000, e.fLnS1000Err);
          if (w) {
            wcdf.push_back(make_pair(s2 / fMaxSin2Theta, w));
            fNEventsAboveCut += w;
          } else if (!fSFunc(e.fLnS1000 - (lnSCut + minLnAtt), e.fLnS1000Err))
            break;
        }
      }
      SortByFirst(wcdf);
      return wcdf;
    }

    template<class Binning>
    std::vector<std::vector<std::pair<double, double>>>
    MakeCDF2D(const Attenuation& att, const Binning& lnS38Binning)
    {
      using namespace std;
      using namespace utl;
      using namespace utl::Accumulator;

      typedef pair<double, double> DD;
      Bin<Collect<DD>, Binning> binLnS38(lnS38Binning);
      const unsigned int nlns = binLnS38.GetNBins();
      const double minLnSCut = binLnS38.GetBinLowerEdge(0);
      const unsigned int ns2 = fSin2ThetaBins.GetNBins();
      const vector<Collect<Event>>& s2Bins = fSin2ThetaBins.GetAllBins();
      for (unsigned int is2 = 0; is2 < ns2; ++is2) {
        const double minLnAtt =
          GetMinLnAttenuation(att, fSin2ThetaBins.GetBinLowerEdge(is2),
                                   fSin2ThetaBins.GetBinUpperEdge(is2));
        for (const Event& e : s2Bins[is2].GetCollection()) {
          const double s2 = e.fSin2Theta;
          const double lnS38 = e.fLnS1000 - log(att(s2));
          const double err = e.fLnS1000Err;
          double wup = fSFunc(lnS38 - binLnS38.GetBinUpperEdge(nlns-1), err);
          for (int ilns = nlns-1; ilns >= 0; --ilns) {
            const double wdn =
              fSFunc(lnS38 - binLnS38.GetBinLowerEdge(ilns), err);
            const double w = wdn - wup;
            if (w)
              binLnS38[ilns](make_pair(s2 / fMaxSin2Theta, w));
            wup = wdn;
          }
          if (!wup && !fSFunc(e.fLnS1000 - (minLnSCut + minLnAtt), err))
            break;
        }
      }

      vector<vector<DD>> ret(nlns);
      for (unsigned int i = 0; i < nlns; ++i) {
        ret[i] = binLnS38[i].GetCollection();
        sort(ret[i].begin(), ret[i].end(),
             [](const DD& a, const DD& b) { return a.first < b.first; });
      }
      ResetCounters();
      for (const auto& row : ret) {
        double rowSum = 0;
        for (const DD& xw : row)
          rowSum += xw.second;
        fNEventsAboveCut += rowSum;
        fLnSMinMaxN(rowSum);
      }

      return ret;
    }

    // Chi^2 measure of uniformity
    double
    GetChi2Ndof(const Attenuation& att, const double lnSCut)
    {
      const std::vector<double> pdf = MakePDF(att, lnSCut);
      return utl::UniformChi2Ndof(pdf, fNEventsAboveCut);
    }

    double
    GetPoissonChi2(const Attenuation& att, const double lnSCut)
    {
      const std::vector<double> pdf = MakePDF(att, lnSCut);
      return utl::UniformPoissonChi2(MakePDF(att, lnSCut), fNEventsAboveCut);
    }

    // Kolmogorov-Smirnov measure of uniformity
    double
    GetKS(const Attenuation& att, const double lnSCut)
    {
      const std::vector<std::pair<double, double>> cdf = MakeCDF(att, lnSCut);
      return utl::UniformKolmogorovSmirnov(cdf, fNEventsAboveCut);
    }

    // Anderson-Drling measure of uniformity
    double
    GetAndersonDarling(const Attenuation& att, const double lnSCut)
    {
      const std::vector<std::pair<double, double>> cdf = MakeCDF(att, lnSCut);
      return utl::UniformAndersonDarling(cdf, fNEventsAboveCut);
    }

    template<class Binning>
    double
    GetPoissonChi22D(const Attenuation& att, const Binning& lnSBinning)
    {
      return utl::UniformPoissonChi2(MakePDF2D(att, lnSBinning));
    }

    template<class Binning>
    double
    GetAndersonDarling2D(const Attenuation& att, const Binning& lnSBinning)
    {
      return utl::UniformAndersonDarling(MakeCDF2D(att, lnSBinning));
    }

    double fMaxSin2Theta = 0;
    utl::Bin<
      BinType,
      utl::EquidistantInNorm<Event>
    > fSin2ThetaBins;
    const SFunction& fSFunc;
    unsigned int fNEvents = 0;
    double fNEventsAboveCut = 0;
    utl::Accumulator::MinMax<double> fSin2ThetaMinMaxN =
      utl::Accumulator::MinMax<double>(std::numeric_limits<double>::infinity(),
                                       -std::numeric_limits<double>::infinity());
    utl::Accumulator::MinMax<double> fLnSMinMaxN =
      utl::Accumulator::MinMax<double>(std::numeric_limits<double>::infinity(),
                                       -std::numeric_limits<double>::infinity());

  };

}


#endif