QuantileCuts.cxx

Go to the documentation of this file.

#include "3FlavorAna/Cuts/QuantileCuts.h"
#include "TH2.h"

#include <iostream>
#include <iomanip>

namespace ana{
  /// \brief Returns a TH2D with xAxis as independentAxis and yAxis as quantileAxis
  /// Quantiles in yAxis can be calculated from this using \ref ana::GetQuantileBins
  TH2* MakeQuantileHistogram(SpectrumLoader& loader,
                             const HistAxis& independentAxis,
                             const HistAxis& quantileAxis,
                             const Cut& cut,
                             const SpillCut& spillCut,
                             const SystShifts& shift,
                             const Var& wei){

    loader.SetSpillCut(spillCut);
    Spectrum spectrum2D(loader, independentAxis, quantileAxis, cut, shift, wei);
    loader.Go();
    return spectrum2D.ToTH2(spectrum2D.POT());

  }

  /// \brief Returns a 2D vector
  /// First index is the independentAxis bin number
  /// Second index is the high bin edge for a quantile
  /// [i][j-1] is jth quantile's low  edge for ith independentAxis bin
  /// [i][j]   is jth quantile's high edge for ith independentAxis bin
  std::vector<std::vector<double>> GetQuantileBins(TH2* quantileHist,
                                                   const HistAxis& independentAxis,
                                                   const HistAxis& quantileAxis,
                                                   const unsigned int& numQuantiles,
                                                   const bool verbose){
    if(independentAxis.NDimensions() != 1 ||
       quantileAxis.NDimensions() != 1){
      std::cout << "Error: GetQuantileBins() not equipped for multi-dimensional axes." << std::endl;
      abort();
    }

    const Binning independentBins = independentAxis.GetBinnings()[0];
    const Binning quantileBins = quantileAxis.GetBinnings()[0];

    if(verbose){
      std::cout << "GetQuantileBins - Calculating Quantile bins" << std::endl;
      std::cout << "independentAxis has " << independentBins.NBins() << " bins and is labelled \'" << independentAxis.GetLabels()[0] << "\'" << std::endl;
      std::cout << "quantileAxis    has " << quantileBins.NBins()    << " bins and is labelled \'" << quantileAxis.GetLabels()[0]    << "\'" << std::endl;
    }

    std::vector<std::vector<double>> xBinQuantilesVec;

    //Copied from Luke Vinton
    //bin=0 is underflow
    //independentBins.NBins()+1 is overflow
    for(int bin = 0; bin <= independentBins.NBins()+1; bin++){

      TH1* hTempProj = quantileHist->ProjectionY(Form("quantileHistProjYBin%d", bin), bin, bin);
      hTempProj->Scale( (double)1 / (double)hTempProj->Integral());
      double probs[numQuantiles];
      double q[numQuantiles];
      for(unsigned int i = 0; i<numQuantiles; i++){
        probs[i] = (i+1) * ((double)1)/ ((double)numQuantiles);
      }

      hTempProj->GetQuantiles(numQuantiles, q, probs);
      //Need to return what the low bin edge is as well, so return a vector of size numQuantiles+1
      std::vector< double > quantilesVec(numQuantiles+1);
      quantilesVec[0] = quantileHist->GetYaxis()->GetBinLowEdge(1);
      for(unsigned int i = 0; i<numQuantiles; i++){
        quantilesVec[i+1] = q[i];
      }

      xBinQuantilesVec.push_back(quantilesVec);
      if(hTempProj) delete hTempProj;
    }//loop over independentAxis bins

    //print out the bin edges
    if(verbose){
      for(unsigned int bin=0;bin<xBinQuantilesVec.size();bin++){
        std::cout << "independentAxis bin: " << std::setw(4) << bin << " "
                  << "lowEdge: " << std::setw(10) << quantileHist->GetXaxis()->GetBinLowEdge(bin) << " "
                  << "centre: "  << std::setw(10) << quantileHist->GetXaxis()->GetBinCenter(bin)  << " "
                  << "quantile bin edges:";
        for(unsigned int edge=0;edge<xBinQuantilesVec[bin].size();edge++) std::cout << " " << std::setw(10) << xBinQuantilesVec[bin][edge];
        std::cout << std::endl;
      }
    }

    return xBinQuantilesVec;
  }

  /// \brief Returns a 2D vector
  /// First index is the independentAxis bin number
  /// Second index is the high bin edge for a quantile
  /// [i][j-1] is jth quantile's low  edge for ith independentAxis bin
  /// [i][j]   is jth quantile's high edge for ith independentAxis bin
  std::vector<std::vector<double>> GetQuantileBins(SpectrumLoader& loader,
                                                   const HistAxis& independentAxis,
                                                   const HistAxis& quantileAxis,
                                                   const Cut& cut,
                                                   const unsigned int& numQuantiles,
                                                   const SpillCut& spillCut,
                                                   const SystShifts& shift,
                                                   const Var& wei,
                                                   const bool verbose){

    TH2* quantileHist = MakeQuantileHistogram(loader, independentAxis, quantileAxis, cut, spillCut, shift, wei);

    std::vector<std::vector<double>> xBinQuantilesVec = GetQuantileBins(quantileHist, independentAxis, quantileAxis, numQuantiles, verbose);

    //clean up after ourselves
    if(quantileHist) delete quantileHist;

    return xBinQuantilesVec;
  }

  /// \brief Class to Generate a Quantile \ref ana::Cut that is specialised via the constructor arguments
  /// \param independentAxis - each bin on this axis will be split into the quantiles (i.e. a range of quantileAxis bins)
  /// \param quantileAxis - the axis to split into quantiles. This is done separately for each independentAxis bin
  /// \param quantile - the quantile to pass, other quantiles fail
  /// \param quantileBins - 2D array specifying quantileAxis quantile edges. 1st index is over independentAxis, second over quantile

  QuantileCutGenerator::QuantileCutGenerator(const HistAxis& independentAxis,
                                             const HistAxis& quantileAxis,
                                             const unsigned int& quantile,
                                             const std::vector<std::vector<double>>& quantileBins):
    fIndependentAxis(independentAxis), fQuantileAxis(quantileAxis), fQuantile(quantile), fQuantileBins(quantileBins)
    {
      if(independentAxis.NDimensions() != 1 ||
         quantileAxis.NDimensions() != 1){
        std::cout << "Error: QuantileCutGenerator not equipped for multi-dimensional axes." << std::endl;
        abort();
      }

      assert(fQuantileBins.size() > 0 && "fQuantileBins must have more than 0 bins in the independent variable");
      //fQuantileBins[] has the TH1 binning schema
      //fQuantileBins[].size = number of fQuantiles + 1 as we keep both the low edge and high edge
      assert((fQuantile > 0 && fQuantile < fQuantileBins[0].size()) && "fQuantile must be > 0 and <= number of fQuantiles");
    }

  bool QuantileCutGenerator::operator() (const caf::SRProxy* sr) const{
    double varX = fIndependentAxis.GetVars()[0](sr);
    double varY = fQuantileAxis.GetVars()[0](sr);
    unsigned int binX = (unsigned int)(fIndependentAxis.GetBinnings()[0].FindBin(varX));

    //fQuantileBins first index is independentAxis bin number, where 0 is underflow
    assert(binX < fQuantileBins.size() && "binX >= fQuantileBins.size()");
    std::vector<double> thisQuantilesBins = fQuantileBins[binX];
    assert(fQuantile > 0 && fQuantile < thisQuantilesBins.size() && "fQuantile <= 0 || fQuantile >= thisQuantilesBins.size()");
    double lowEdge=thisQuantilesBins[fQuantile-1];
    double highEdge=thisQuantilesBins[fQuantile];
    return varY >= lowEdge && varY < highEdge;
  }


  /// \brief Returns a Cut that will only pass events falling into a particular quantile of quantileAxis, where quantiles have been calculated for each independentAxis bin
  /// Quantiles have been precalculated in 'quantileBins' (number of quantiles, bin edges etc.)
  const Cut QuantileCut(const std::vector<std::vector<double>>& quantileBins,
                        const HistAxis& independentAxis,
                        const HistAxis& quantileAxis,
                        const unsigned int& quantile){
    if(independentAxis.NDimensions() != 1 ||
       quantileAxis.NDimensions() != 1){
      std::cout << "Error: QuantileCut() not equipped for multi-dimensional axes." << std::endl;
      abort();
    }

    assert((quantileBins.size() > 0) && "quantileBins.size() must be > 0");
    unsigned int numQuantiles = quantileBins[0].size();
    assert((quantile > 0 && quantile <= numQuantiles) && "quantile must be > 0 and <= number of quantiles");

    return Cut(QuantileCutGenerator(independentAxis, quantileAxis, quantile, quantileBins));

  }

  /// \brief Returns a vector of Cuts, each one for a different quantile
  /// An individual cut will only pass events falling into a particular quantile of quantileAxis, where quantiles have been calculated for each independentAxis bin
  std::vector<Cut> QuantileCuts(SpectrumLoader& loader,
                                const HistAxis& independentAxis,
                                const HistAxis& quantileAxis,
                                const Cut& cut,
                                const unsigned int& numQuantiles,
                                const SpillCut& spillCut,
                                const SystShifts& shift,
                                const Var& wei,
                                const bool verbose){

    std::vector<std::vector<double>> quantileBins = GetQuantileBins(loader, independentAxis, quantileAxis, cut, numQuantiles, spillCut, shift, wei, verbose);

    std::vector<Cut> cutVec;
    for(unsigned int quantile=1;quantile<=numQuantiles;quantile++){
      cutVec.push_back(QuantileCut(quantileBins, independentAxis, quantileAxis, quantile));
    }
    return cutVec;
  }

  /// \breif: Do the same as the QuantileCuts function but taking in the TH2 instead of making it.
  std::vector<Cut> QuantileCutsFromTH2(TH2* quantileHist,
                                       const HistAxis& independentAxis,
                                       const HistAxis& quantileAxis,
                                       const unsigned int& numQuantiles,
                                       const bool verbose){

    std::vector<std::vector<double>> quantileBins = GetQuantileBins(quantileHist, independentAxis, quantileAxis, numQuantiles, verbose);

    std::vector<Cut> cutVec;
    for(unsigned int quantile=1;quantile<=numQuantiles;quantile++){
      cutVec.push_back(QuantileCut(quantileBins, independentAxis, quantileAxis, quantile));
    }
    return cutVec;
  }

  std::vector<Cut> QuantileAndPIDCutsFromTH2(TH2* quantileHist,
                                       const HistAxis& independentAxis,
                                       const HistAxis& quantileAxis,
                                       const unsigned int& numQuantiles,
               const bool rhc,
                                       const bool verbose){

    std::vector<std::vector<double>> quantileBins = GetQuantileBins(quantileHist, independentAxis, quantileAxis, numQuantiles, verbose);

    std::vector<Cut> cutVec;
    for(unsigned int quantile=1;quantile<=numQuantiles;quantile++){
      cutVec.push_back(QuantileCut(quantileBins, independentAxis, quantileAxis, quantile));
    }
    if (cutVec.size()==4 && rhc){
      cutVec[0] = cutVec[0] && k2018PIDs(0.47, 0.40, 0.80);
      cutVec[1] = cutVec[1] && k2018PIDs(0.48, 0.40, 0.65);
      cutVec[2] = cutVec[2] && k2018PIDs(0.56, 0.40, 0.60);
      cutVec[3] = cutVec[3] && k2018PIDs(0.59, 0.65, 0.50);
    }
    else if (cutVec.size()==4 && !rhc){
      cutVec[0] = cutVec[0] && k2018PIDs(0.56, 0.60, 0.55);
      cutVec[1] = cutVec[1] && k2018PIDs(0.56, 0.60, 0.55);
      cutVec[2] = cutVec[2] && k2018PIDs(0.58, 0.60, 0.55);
      cutVec[3] = cutVec[3] && k2018PIDs(0.59, 0.60, 0.55);
    }
    else for(unsigned int quant=0;quant<numQuantiles;quant++) cutVec[quant] = cutVec[quant] && k2018PIDs();

    return cutVec;
  }




}//namespace