CovarianceBinUtility.cxx

Go to the documentation of this file.

//
// Created by Brian Rebel on 11/29/18.
//

#include "messagefacility/MessageLogger/MessageLogger.h"
#include "cetlib_except/exception.h"

#include "NovaDAQConventions/DAQConventions.h"
#include "CovarianceMatrixFit/utilities/CovarianceBinUtility.h"
#include "CovarianceMatrixFit/dataProducts/Structs.h"
#include "CovarianceMatrixFit/dataProducts/StaticFuncs.h"

namespace cmf {

  static CovarianceBinUtility* gCovBinUtil = nullptr;

  //-----------------------------------------------------------------------------
  CovarianceBinUtility* CovarianceBinUtility::Instance()
  {
    if(gCovBinUtil == nullptr) gCovBinUtil = new CovarianceBinUtility();

    return gCovBinUtil;
  }
  
  //-----------------------------------------------------------------------------
  CovarianceBinUtility::CovarianceBinUtility()
  {

  }

  //-----------------------------------------------------------------------------
  int CovarianceBinUtility::NumSelectionBins(cmf::SelectionType_t const& selType,
                                             cmf::DetType_t       const& detType)
  {
    if(cmf::IsNuESelected(selType)){
      if(selType == cmf::kNuESelectionPeripheral) return fNuEPeriphHighEdges.size();
      return (detType == cmf::kNEARDET) ? fNuENDHighEdges.size() : fNuEFDHighEdges.size();
    }
    else if(cmf::IsNuMuSelected(selType)){
      return fNuMuHighEdges.size();
    }
    else if(selType == cmf::kNCSelection){
      return (detType == cmf::kNEARDET) ? fNCNDHighEdges.size() : fNCFDHighEdges.size();
    }

    throw cet::exception("CovarianceBinUtility")
        << "selection type "
        << selType
        << " is not recognized";

    return 0;
  }

  //-----------------------------------------------------------------------------
  void CovarianceBinUtility::FillOffsetAndBinMaps(std::map<long, int>                       & offsetMap,
                                                  std::vector<long>                         & binMap,
                                                  cmf::SelectionUtility::DetBeamSelSet const& detBeamSelSet)
   {
     long key;
     int  offset = 0;
     int  bins;

     binMap   .clear();
     offsetMap.clear();
     
     // loop over detectors, beam type and selection to fill the maps
     for(auto const& dbsItr : detBeamSelSet){
       for(auto const& selItr : dbsItr.Selections() ){

         // there is no MidPID selection for nues any more and no nue peripheral
         // selection in the nd
         if(dbsItr.Detector() == cmf::kNEARDET &&
            selItr == cmf::kNuESelectionPeripheral)
           continue;

         // set the file type to 0 so that it does not contribute to the
         // key value in this case
         key = DetectorBeamSelectionTypesToKey(dbsItr.Detector(),
                                               dbsItr.BeamType(),
                                               selItr);

         offsetMap[key] = offset;

         bins = this->NumSelectionBins(selItr, dbsItr.Detector());
         for(int b = offset; b < bins + offset; ++b){
           binMap.emplace_back(key);
           LOG_DEBUG("CovarianceBinUtility")
               << cmf::KeyToString(key)
               << " key "
               << key
               << " offset "
               << offsetMap.find(key)->second
               << " "
               << b
               << " "
               << binMap[b];
         }

         offset += bins;
       } // end loop over selection types
     } // end loop over input DetBeamSels

   }

  //---------------------------------------------------------------------------
  // create map from the vectors of bin edges read from fhicl
  std::map<double, int> CovarianceBinUtility::EnergyBinMap(std::vector<double> const& vec)
  {
    std::map<double, int> ebMap;
    for (size_t i = 0; i < vec.size(); ++i)
      ebMap[vec.at(i)] = i;

    return ebMap;
  }

  //---------------------------------------------------------------------------
  // helper method to setup the energy-to-bin maps 
  void CovarianceBinUtility::SetBinning(fhicl::ParameterSet const& fBinConfig, cmf::SelectionType_t const& selType)
  {
    auto ndConfig = fBinConfig.get<fhicl::ParameterSet>("nd");
    auto fdConfig = fBinConfig.get<fhicl::ParameterSet>("fd");

    // numu selection has the same binning in nd and fd
    // TODO if the numu selection is ever not the same in the detectors we
    // need to update this
    if(cmf::IsNuMuSelected(selType)){
      fNuMuLowEdge   = ndConfig.get<double>("low_edge");
      fNuMuHighEdges = this->EnergyBinMap(ndConfig.get<std::vector<double>>("high_edges"));
    }
    else if (cmf::IsNuESelected(selType)){
      fNuENDLowEdge   = ndConfig.get<double>("low_edge");
      fNuENDHighEdges = this->EnergyBinMap(ndConfig.get<std::vector<double>>("high_edges"));
      fNuEFDLowEdge   = fdConfig.get<double>("low_edge");
      fNuEFDHighEdges = this->EnergyBinMap(fdConfig.get<std::vector<double>>("high_edges"));

      // nue selection also has peripheral sample at the FD
      auto peConfig       = fBinConfig.get<fhicl::ParameterSet>("peripheral");
      fNuEPeriphLowEdge   = peConfig.get<double>("low_edge");
      fNuEPeriphHighEdges = this->EnergyBinMap(peConfig.get<std::vector<double>>("high_edges"));
    }
    else {
      fNCNDLowEdge   = ndConfig.get<double>("low_edge");
      fNCNDHighEdges = this->EnergyBinMap(ndConfig.get<std::vector<double>>("high_edges"));
      fNCFDLowEdge   = fdConfig.get<double>("low_edge");
      fNCFDHighEdges = this->EnergyBinMap(fdConfig.get<std::vector<double>>("high_edges"));
    }

  }

  //-----------------------------------------------------------------------------
  void CovarianceBinUtility::Initialize(fhicl::ParameterSet const& pset)
  {

    // grab configurations from parameter set
    auto fNuMuEdgeConfig = pset.get< fhicl::ParameterSet >("numu_bin_edges");
    auto fNuEEdgeConfig  = pset.get< fhicl::ParameterSet >("nue_bin_edges");
    auto fNCEdgeConfig   = pset.get< fhicl::ParameterSet >("nc_bin_edges");

    // using the pre-existing selection enum, doesn't matter
    // since we just need to differentiate between numu nue and nc
    this->SetBinning(fNuMuEdgeConfig, cmf::kNuMuSelection);
    this->SetBinning(fNuEEdgeConfig , cmf::kNuESelection);
    this->SetBinning(fNCEdgeConfig  , cmf::kNCSelection);

    LOG_DEBUG("CovarianceBinUtility")
      << "\nnum NuMuBins:     "  << fNuMuHighEdges.size()      
      << "\nnum NuENDBins:    "  << fNuENDHighEdges.size()     
      << "\nnum NuEFDBins:    "  << fNuEFDHighEdges.size()     
      << "\nnum NuEPeriphBins:"  << fNuEPeriphHighEdges.size() 
      << "\nnum NCNDBins:     "  << fNCNDHighEdges.size()      
      << "\nnum NCFDBins:     "  << fNCFDHighEdges.size();

    FillBinToEnergyMap(fNuMuHighEdges,
                       fNuMuBinToEnergy,
                       fNuMuLowEdge);

    FillBinToEnergyMap(fNuENDHighEdges,
                       fNuENDBinToEnergy,
                       fNuENDLowEdge);

    FillBinToEnergyMap(fNuEFDHighEdges,
                       fNuEFDBinToEnergy,
                       fNuEFDLowEdge);

    FillBinToEnergyMap(fNuEPeriphHighEdges,
                       fNuEPeriphBinToEnergy,
                       fNuEPeriphLowEdge);

    FillBinToEnergyMap(fNCNDHighEdges,
                       fNCNDBinToEnergy,
                       fNCNDLowEdge);

    FillBinToEnergyMap(fNCFDHighEdges,
                       fNCFDBinToEnergy,
                       fNCFDLowEdge);

    // call the function to check that the mapping from 
    // energy to bin and back again is consistent
    //this->CheckBinToEnergyMapping();

    LOG_DEBUG("CovarianceBinUtility")
    << "Fill All Selection Offset and Bin maps";

    this->FillOffsetAndBinMaps(fAllSelsOffsetMap,
                               fAllSelsBinToKey,
                               cmf::SelectionUtility::Instance()->AllAvailSelections());

    LOG_DEBUG("CovarianceBinUtility")
        << "Fill run time Selection Offset and Bin maps";

    this->FillOffsetAndBinMaps(fOffsetMap,
                               fBinToKey,
                               cmf::SelectionUtility::Instance()->SelectionsToUse());

    LOG_VERBATIM("CovarianceBinUtility")
      << "Using "
      << fBinToKey.size()
      << " bins for analysis out of a possible "
      << fAllSelsBinToKey.size();
  }

  //-----------------------------------------------------------------------------
  CovarianceBinUtility::~CovarianceBinUtility()
  {
    fNuMuHighEdges   .clear();
    fNuENDHighEdges  .clear();
    fNuEFDHighEdges  .clear();
    fNCHighEdges     .clear();
    fNCNDHighEdges   .clear();
    fNCFDHighEdges   .clear();
    fAllSelsOffsetMap.clear();
    fOffsetMap       .clear();
    fBinToKey        .clear();
    fAllSelsBinToKey .clear();

    fNuEPeriphHighEdges  .clear();
    fNuMuBinToEnergy     .clear();
    fNuENDBinToEnergy    .clear();
    fNuEFDBinToEnergy    .clear();
    fNuEPeriphBinToEnergy.clear();
    fNCNDBinToEnergy     .clear();
    fNCFDBinToEnergy     .clear();
 
  }

  //-----------------------------------------------------------------------------
  void CovarianceBinUtility::FillBinToEnergyMap(std::map<double, int> const& selectionHighEdges,
                                                std::map<int, double>      & selectionBinToEnergy,
                                                double                       lowestBound)
  {
    // fill the bin to energy maps for fast look up
    double lowerBound = lowestBound;
    for(auto const& itr : selectionHighEdges){
      selectionBinToEnergy[itr.second] = lowerBound + (itr.first - lowerBound) * 0.5;
      
      LOG_DEBUG("CovarianceBinUtility")
          << "FillBinToEnergyMap energy "
          << itr.first
          << " bin "
          << itr.second
          << " bin to energy "
          << selectionBinToEnergy.find(itr.second)->second;

      lowerBound = itr.first;
    }
  }
  
  //-----------------------------------------------------------------------------
//  void CovarianceBinUtility::CheckBinToEnergyMapping()
//  {
//    // loop over the selections and detectors to get the correct number of
//    // bins for each and then get the energy for each bin.  Use that value to
//    // find the bin again and compare to the loop iterator
//    std::set<cmf::SelectionType_t> selectionTypes({kNuESelectionLowPID,
//                                                   kNuESelectionHighPID,
//                                                   kNuESelectionPeripheral,
//                                                   kNuMuSelectionQ1,
//                                                   kNuMuSelectionQ2,
//                                                   kNuMuSelectionQ3,
//                                                   kNuMuSelectionQ4,
//                                                   kNCSelection});
//    std::set<cmf::DetType_t> detectors({cmf::kNEARDET,
//                                             cmf::kFARDET});
//    double binEnergy;
//    size_t binNumber;
//
//    // create a MetaData for making the following easier,
//    // set the epoch to be 1, ie  epoch 1, FHC
//    cmf::MetaData md;
//    md.period = 1;
//
//    for(auto const& st : selectionTypes){
//      md.selectionType = st;
//
//      for(auto const& det : detectors){
//        md.detector = det;
//
//        auto const& highEdges = this->SelectionHighEdges(md);
//        for(auto const& itr : highEdges){
//
//          // no peripheral selection in the ND
//          if(det == cmf::kNEARDET &&
//             st  == cmf::kNuESelectionPeripheral) continue;
//
//          binEnergy = this->BinToEnergy(itr.second + this->KeyToOffset(md.DetectorBeamSelectionKey()), md);
//          binNumber = this->EnergyToBin(binEnergy, md);
//
//          LOG_DEBUG("CovarianceBinUtility")
//              << cmf::cSelectionType_Strings[st]
//              << " "
//              << md.DetectorString()
//              << " look up of bin number:  "
//              << itr.second
//              << " found bin number "
//              << binNumber
//              << " offset "
//              << this->KeyToOffset(md.DetectorBeamSelectionKey())
//              << " found energy "
//              << binEnergy;
//
//        }// end loop over bins
//      } // end loop over detectors
//    } // end loop over selections
//  }
  
  //......................................................................
  // There are 4 quantiles for numu events, 2 PID sets + 1 peripheral bin
  // for nue events, and 1 set for NC events.  Multiply that by 2 for the
  // number of horn configurations
  size_t CovarianceBinUtility::TotalBins(bool allSels)
  {

    LOG_DEBUG("CovarianceBinUtility")
      << " there are numu (each quartile): "
      << fNuMuHighEdges.size()
      << " nue low: "
      << fNuENDHighEdges.size()
      << " nue high: "
      << fNuEFDHighEdges.size()
      << " nue peripheral: "
      << fNuEPeriphHighEdges.size()
      << " nc ND: "
      << fNCNDHighEdges.size()
      << " nc FD: "
      << fNCFDHighEdges.size()
      << " total bins for all selections: "
      << fAllSelsBinToKey.size()
      << " total bins for run time selections: "
      << fBinToKey.size()
      << " "
      << allSels;
    
    if(allSels) return fAllSelsBinToKey.size();

    return fBinToKey.size();    
  }
 
  //......................................................................
  double const& CovarianceBinUtility::SelectionLowEdge(cmf::MetaData const& md)
  {
    if(md.IsNuMuSelected())
      return this->NuMuLowEdge();
    else if(md.IsNuESelected()){
      return this->NuELowEdge(md);
    }
    
    return this->NCLowEdge(md.detector);
  }

  //......................................................................
  double const& CovarianceBinUtility::NuELowEdge(cmf::MetaData const& md) 
  {
    if(md.detector == cmf::kNEARDET) return fNuENDLowEdge;
    if(md.selectionType == cmf::kNuESelectionPeripheral) return fNuEPeriphLowEdge;

    return fNuEFDLowEdge;
  }

  //......................................................................
  double const& CovarianceBinUtility::NCLowEdge(cmf::DetType_t const& det)
  {
    if(det == cmf::kNEARDET) return fNCNDLowEdge;
    
    return fNCFDLowEdge;
  }

  //......................................................................
  std::map<double, int> const& CovarianceBinUtility::NuEHighEdges(cmf::DetType_t const& det)
  {
    if(det == cmf::kNEARDET) return fNuENDHighEdges;
    
    return fNuEFDHighEdges;
  }

  //......................................................................
  std::map<double, int> const& CovarianceBinUtility::NCHighEdges(cmf::DetType_t const& det)
  {
    if(det == cmf::kNEARDET) return fNCNDHighEdges;
    
    return fNCFDHighEdges;
  }

  //......................................................................
  std::map<double, int> const& CovarianceBinUtility::SelectionHighEdges(cmf::MetaData const& md)
  {
    if(md.IsNuMuSelected())
      return this->NuMuHighEdges();
    else if(md.IsNuESelected()){
      if(md.selectionType == cmf::kNuESelectionPeripheral) return fNuEPeriphHighEdges;
      return this->NuEHighEdges(md.detector);
    }
    
    return this->NCHighEdges(md.detector);
  }

  //......................................................................
  std::map<double, int> const& CovarianceBinUtility::SelectionHighEdges(long const& key)
  {
    auto sel = cmf::KeyToSelectionType(key);
    auto det = cmf::KeyToDetectorType(key);

    if(cmf::IsNuMuSelected(sel))
      return this->NuMuHighEdges();
    else if(cmf::IsNuESelected(sel)){
      if(sel == cmf::kNuESelectionPeripheral) return fNuEPeriphHighEdges;
      return this->NuEHighEdges(det);
    }

    return this->NCHighEdges(det);
  }

  //......................................................................
  void CovarianceBinUtility::SelectionHistBinning(cmf::MetaData  const& md,
                                                  std::vector<double> & bins)
  {
    std::vector<double> tmpbins;
    auto const& highEdgeMap = this->SelectionHighEdges(md);
    tmpbins.push_back(this->SelectionLowEdge(md));

    for(auto const& itr : highEdgeMap) tmpbins.push_back(itr.first);

    std::swap(tmpbins, bins);
  }

  //......................................................................
  int CovarianceBinUtility::KeyToOffset(long const& key,
                                        bool        allSels)
  {
    if(allSels) return this->KeyToOffset(key, fAllSelsOffsetMap);
    return this->KeyToOffset(key, fOffsetMap);
  }

  //......................................................................
  int CovarianceBinUtility::KeyToOffset(long                const& key,
                                        std::map<long, int> const& offsetMap)
  {

    if(offsetMap.count(key) < 1)
      throw cet::exception("CovarianceBinUtility")
          << "could not find offset in map for key "
          << key
          << " " 
          << cmf::cBeamType_Strings[cmf::KeyToBeamType(key)]
          << " "
          << cmf::cDetType_Strings[cmf::KeyToDetectorType(key)]
          << " " 
          << cmf::cSelectionType_Strings[cmf::KeyToSelectionType(key)];

    return offsetMap.find(key)->second;
  }

  //......................................................................
  int CovarianceBinUtility::EnergyToHistogramBin(double        const& energy,
                                                 cmf::MetaData const& md)
  {
    return this->EnergyToBin(energy, md) + 1;
  }

  //......................................................................
  int CovarianceBinUtility::EnergyToBin(double        const& energy,
                                        cmf::MetaData const& md)
  {
    auto const& edgeMap = this->SelectionHighEdges(md);
    
    auto const& itr = edgeMap.upper_bound(energy);
    if(itr != edgeMap.end()) return itr->second;
    
    return 0;
  }

  //......................................................................
  std::map<int, double> const& CovarianceBinUtility::SelectionBinToEnergy(long key)
  {
    auto sel = cmf::KeyToSelectionType(key);
    auto det = cmf::KeyToDetectorType(key);

    if(sel == cmf::kUnknownSelection ||
       det == cmf::kUnknownDet )
      throw cet::exception("Unknown Key")
      << "key "
      << key
      << "corresponds to "
      << cmf::cSelectionType_Strings[sel]
      << " and "
      << cmf::cDetType_Strings[det];

    if(cmf::IsNuESelected(sel)){
      if(sel == cmf::kNuESelectionPeripheral)  return fNuEPeriphBinToEnergy;
      if(det == cmf::kNEARDET) return fNuENDBinToEnergy;
      return fNuEFDBinToEnergy;
    }
    else if(cmf::IsNuMuSelected(sel)){
      return fNuMuBinToEnergy;
    }
    else{
      if (det == cmf::kNEARDET) return fNCNDBinToEnergy;
      else return fNCFDBinToEnergy;
    }
  }

  //-----------------------------------------------------------------------------
  cmf::DetType_t CovarianceBinUtility::BinToDetector(int const& bin,
                                                          bool       allSels)
  {
    return cmf::KeyToDetectorType(this->BinToKey(bin, allSels));
  }

  //-----------------------------------------------------------------------------
  cmf::SelectionType_t CovarianceBinUtility::BinToSelection(int const& bin,
                                                            bool       allSels)
  {
    return cmf::KeyToSelectionType(this->BinToKey(bin, allSels));
  }

  //-----------------------------------------------------------------------------
  cmf::BeamType_t CovarianceBinUtility::BinToBeamType(int const& bin,
                                                      bool       allSels)
  {
    return cmf::KeyToBeamType(this->BinToKey(bin, allSels));
  }

  //-----------------------------------------------------------------------------
  long CovarianceBinUtility::BinToKey(int const& bin,
                                      bool       allSels)
  {
    if(allSels) return fAllSelsBinToKey[bin];
    return fBinToKey[bin];
  }

  //......................................................................
  double CovarianceBinUtility::BinToEnergy(int  const& bin,
                                           bool        allSels)
  {
    auto key    = this->BinToKey   (bin, allSels);
    int  offset = this->KeyToOffset(key, allSels);

    auto const& binToEnergy = this->SelectionBinToEnergy(key);

    auto const btoEItr = binToEnergy.find(bin - offset);

    if(btoEItr == binToEnergy.end())
      throw cet::exception("CovarianceBinUtility")
          << "bin "
          << bin
          << " offset "
          << offset
          << " "
          << bin - offset
          << " is not in bin to energy map for key "
          << key
          << " "
          << cmf::KeyToString(key);

    LOG_DEBUG("CovarianceBinUtility")
        << "BinToEnergy: "
        << cmf::KeyToString(key)
        << " bin is "
        << bin
        << " offset "
        << offset
        << " energy "
        << btoEItr->second;

    return btoEItr->second;
  }

  //......................................................................
  std::vector<int> CovarianceBinUtility::EmptyBins(TMatrixD mat){

    std::vector<int> emptyBins;

    for (int i = 0; i < mat.GetNrows(); i++){

      // for every column, check every row is above some
      // tolerance. If it's not, count it as an empty bin
      double rowsum = 0;
      for (int j = 0; j < mat.GetNcols(); j++)
        rowsum += std::abs(mat(i,j)); 

      if (rowsum == 0)
        emptyBins.push_back(i);
    }

    return emptyBins;
  }

  //......................................................................
  TMatrixD CovarianceBinUtility::CompactMatrix(TMatrixD mat, std::vector<int> emptyBins){

    int nRowsAndColumns = mat.GetNcols() - emptyBins.size();
    TMatrixD compactMatrix(nRowsAndColumns, nRowsAndColumns);
    TMatrixD semiCompactMatrix(nRowsAndColumns, mat.GetNcols());

    // correct in i
    int itest = 0;
    for (int i = 0; i < nRowsAndColumns; ++i){
      
      bool isXEmptyBin = false;
      for (size_t ieb = 0; ieb < emptyBins.size(); ++ieb){
        if (emptyBins.at(ieb) == itest)
          isXEmptyBin = true;
      }

      for (int j = 0; j < mat.GetNcols(); ++j){
        semiCompactMatrix(i,j) = mat(itest, j);
      }

      if (isXEmptyBin)
        --i;

      ++itest;

    }

    // correct in i
    int jtest = 0;
    for (int j = 0; j < nRowsAndColumns; ++j){
      
      bool isYEmptyBin = false;
      for (size_t ieb = 0; ieb < emptyBins.size(); ++ieb){
        if (emptyBins.at(ieb) == jtest)
          isYEmptyBin = true;
      }

      for (int i = 0; i < semiCompactMatrix.GetNrows(); ++i){
        compactMatrix(i,j) = semiCompactMatrix(i, jtest);
      }

      if (isYEmptyBin)
        --j;

      ++jtest;

    }

    return compactMatrix;

  }

} // cmf namespace