PlotUtilities.cxx

Go to the documentation of this file.

//
// Created by Brian Rebel on 11/19/19.
//
#include <sstream>
#include <iomanip>

#include "TCanvas.h"
#include "TMath.h"
#include "TLegend.h"
#include "TMarker.h"
#include "TList.h"
#include "TLine.h"
#include "TROOT.h"
#include "TKey.h"
#include "TObjArray.h"

#include "CovarianceMatrixFit/utilities/PlotUtilities.h"
#include "CovarianceMatrixFit/utilities/CovarianceBinUtility.h"
#include "CovarianceMatrixFit/dataProducts/StaticFuncs.h"

namespace cmf{

  static PlotUtilities* gPlotUtil = nullptr;

  //-----------------------------------------------------------------------------
  PlotUtilities* PlotUtilities::Instance()
  {
    if(gPlotUtil == nullptr) gPlotUtil = new PlotUtilities();

    return gPlotUtil;
  }

  //-----------------------------------------------------------------------------
  PlotUtilities::PlotUtilities()
  : fXParamEnum        (cmf::kNumOscParams)
  , fYParamEnum        (cmf::kNumOscParams)
  , fFreeParamEnum     (cmf::kNumOscParams)
  , fXParamExtrema     (std::make_pair(cmf::kGarbageDouble, std::numeric_limits<double>::lowest()))
  , fYParamExtrema     (std::make_pair(cmf::kGarbageDouble, std::numeric_limits<double>::lowest()))
  , fXParamScale       (std::make_pair(1., ""))
  , fYParamScale       (std::make_pair(1., ""))
  , fContourSigmaLabels(std::vector<std::string>(3))
  {
  }

  //-----------------------------------------------------------------------------
  void PlotUtilities::Initialize(fhicl::ParameterSet const& plottingPars)
  {
    fUseTrig                = plottingPars.get<bool       >("UseTrig",               false   );
    fContourLevel1Sigma1D   = plottingPars.get<double     >("ContourLevel1Sigma1D",  1       );
    fContourLevel2Sigma1D   = plottingPars.get<double     >("ContourLevel2Sigma1D",  4       );
    fContourLevel3Sigma1D   = plottingPars.get<double     >("ContourLevel3Sigma1D",  9       );
    fContourLevel1Sigma2D   = plottingPars.get<double     >("ContourLevel1Sigma2D",  2.3     );
    fContourLevel2Sigma2D   = plottingPars.get<double     >("ContourLevel2Sigma2D",  6.18    );
    fContourLevel3Sigma2D   = plottingPars.get<double     >("ContourLevel3Sigma2D",  11.83   );
    fContourSigmaLabels[0]  = plottingPars.get<std::string>("Contour1SigmaLabel",   "1#sigma");
    fContourSigmaLabels[1]  = plottingPars.get<std::string>("Contour2SigmaLabel",   "2#sigma");
    fContourSigmaLabels[2]  = plottingPars.get<std::string>("Contour3SigmaLabel",   "3#sigma");
    fNormalizeBins          = plottingPars.get<bool       >("NormalizeBinWidths",    false   );
    fNuMuNormVal            = plottingPars.get<double     >("NuMuNormVal"       ,    0.1     );
    fNuENormVal             = plottingPars.get<double     >("NuENormVal"        ,    0.5     );
    fNCNormVal              = plottingPars.get<double     >("NCNormVal"         ,    1.0     );
  }

  //-----------------------------------------------------------------------------
  void PlotUtilities::Initialize(fhicl::ParameterSet       const& plottingPars,
                                 cmf::OscParm_t            const& xParEnum,
                                 cmf::OscParm_t            const& yParEnum,
                                 size_t                    const& numXParPoints,
                                 size_t                    const& numYParPoints,
                                 std::pair<double, double> const& parXExtrema,
                                 std::pair<double, double> const& parYExtrema,
                                 bool                      const& logXAxis,
                                 bool                      const& logYAxis)
  {
    this->Initialize(plottingPars);

    fXParamExtrema = parXExtrema;
    fYParamExtrema = parYExtrema;
    fXParamDivs    = numXParPoints;
    fYParamDivs    = numYParPoints;
    fLogXAxis      = logXAxis;
    fLogYAxis      = logYAxis;
    fXParamEnum    = xParEnum;
    fYParamEnum    = yParEnum;

    if(fXParamEnum == fYParamEnum)
      throw cet::exception("ContourMaker")
          << "x and y parameters are the same "
          << cmf::cOscParams_Strings[fXParamEnum]
          << " - that cannot be good, bail";

    fXParamScale = std::make_pair(cmf::cOscParams_Scales[fXParamEnum],
                                  cmf::cOscParams_LatexScales[fXParamEnum]);
    fYParamScale = std::make_pair(cmf::cOscParams_Scales[fYParamEnum],
                                  cmf::cOscParams_LatexScales[fYParamEnum]);

    MF_LOG_VERBATIM("PlotUtilities")
        << "Plot Utilities: "
        << "X Parameter: "
        << cmf::cOscParams_Strings[fXParamEnum]
        << " ("
        << fXParamExtrema.first
        << ", "
        << fXParamExtrema.second
        << ") numDivs: "
        << fXParamDivs
        << " Log axis: "
        << fLogXAxis
        << " Y Parameter: "
        << cmf::cOscParams_Strings[fYParamEnum]
        << " ("
        << fYParamExtrema.first
        << ", "
        << fYParamExtrema.second
        << ") numDivs: "
        << fYParamDivs
        << " Log axis: "
        << fLogYAxis;

    float paramDiv;
    float paramStart;
    float paramStop;

    // set the size of the vectors for the parameter bins to be the number of
    // divisions + 2 to keep track of the end points
    fXParamBinEdges = std::vector<float>(fXParamDivs + 2);
    fYParamBinEdges = std::vector<float>(fYParamDivs + 2);

    if(fLogXAxis){
      // convert extrema to log values - if we are logging a parameter we
      // have a scale of 1, so don't bother using the f*ParamScale
      float logXParamExtremaLow  = std::log10(fXParamExtrema.first);
      float logXParamExtremaHigh = std::log10(fXParamExtrema.second);
      paramDiv   = (logXParamExtremaHigh - logXParamExtremaLow)/float(fXParamDivs);
      paramStart = (logXParamExtremaLow  - 0.5 * paramDiv);
      paramStop  = (logXParamExtremaHigh + 0.5 * paramDiv);
      for(size_t i = 0; i < fXParamBinEdges.size() - 1; ++i){
        fXParamBinEdges[i] = std::pow(10,(paramStart + (i * paramDiv)));
      }
      fXParamBinEdges.back() = std::pow(10, paramStop);
    }
    else {

      paramDiv   = (fXParamExtrema.second - fXParamExtrema.first )/float(fXParamDivs);
      paramStart = (fXParamExtrema.first  - 0.5 * paramDiv) * fXParamScale.first;
      paramStop  = (fXParamExtrema.second + 0.5 * paramDiv) * fXParamScale.first;
      for (size_t i = 0; i < fXParamBinEdges.size() - 1; ++i){
        fXParamBinEdges[i] = paramStart + (i * paramDiv * fXParamScale.first);
      }
      fXParamBinEdges.back() = paramStop;
    }

    // for(auto const& itr : fXParamBinEdges){
    //   MF_LOG_VERBATIM("PlotUtilities")
    //  << "x bin edges: "
    //  << itr;
    // }

    if(fLogYAxis)
    {
      // convert extrema to log values - if we are logging a parameter we
      // have a scale of 1, so don't bother using the f*ParamScale
      float logYParamExtremaLow  = std::log10(fYParamExtrema.first);
      float logYParamExtremaHigh = std::log10(fYParamExtrema.second);
      paramDiv   = (logYParamExtremaHigh - logYParamExtremaLow)/float(fYParamDivs);
      paramStart = (logYParamExtremaLow  - 0.5 * paramDiv);
      paramStop  = (logYParamExtremaHigh + 0.5 * paramDiv);
      for (size_t i = 0; i < fYParamBinEdges.size() - 1; ++i){
        fYParamBinEdges[i] = std::pow(10,(paramStart + (i * paramDiv)));
      }
      fYParamBinEdges.back() = std::pow(10, paramStop);
    }
    else {
      paramDiv   = (fYParamExtrema.second - fYParamExtrema.first )/float(fYParamDivs);
      paramStart = (fYParamExtrema.first  - 0.5 * paramDiv) * fYParamScale.first;
      paramStop  = (fYParamExtrema.second + 0.5 * paramDiv) * fYParamScale.first;

      for (size_t i = 0; i < fYParamBinEdges.size() - 1; ++i){
        fYParamBinEdges[i] = paramStart + (i * paramDiv * fYParamScale.first);
      }
      fYParamBinEdges.back() = paramStop;
    }

    // for(auto const& itr : fYParamBinEdges){
    //   MF_LOG_VERBATIM("PlotUtilities")
    //  << "y bin edges: "
    //  << itr;
    // }

  }

  //......................................................................
  void PlotUtilities::Make1DPlot(art::TFileDirectory                   & tfd,
                                 std::string                      const& chiSqrName,
                                 std::map<double, double>         const& paramChiSqr,
                                 bool                                    xPar)
  {

    std::vector<float>             paramBins  = fYParamBinEdges;
    cmf::OscParm_t                 paramEnum  = fYParamEnum;
    std::pair<double, std::string> paramScale = fYParamScale;
    if(xPar){
      paramBins  = fXParamBinEdges;
      paramEnum  = fXParamEnum;
      paramScale = fXParamScale;
    }

    std::string grName = chiSqrName + cmf::cOscParams_Strings[paramEnum] + "DeltaChiSqr";

    //make vectors holding the points for the graph
    std::vector<double> param;
    std::vector<double> deltaChiSqr;

    for(auto itr : paramChiSqr){
      param      .emplace_back(itr.first * paramScale.first);
      deltaChiSqr.emplace_back(itr.second);

      MF_LOG_DEBUG("PlotUtilities")
      << grName
      << " "
      << param.back()
      << " "
      << deltaChiSqr.back();
    }

    // make the graph
    auto *gr = tfd.makeAndRegister<TGraph>(grName.c_str(),
                                           grName.c_str(),
                                           param.size(),
                                           param.data(),
                                           deltaChiSqr.data());

    // make a histogram to set the axis ranges, maximum DeltaChiSqr = 15
    std::string histName = grName + "Hist";
    std::string title    = ";" + cmf::cOscParams_Strings_Latex[paramEnum];
    if     (paramEnum == cmf::kTh23 && fUseTrig) title = ";sin^{2}(#theta_{23})";
    else if(paramEnum == cmf::kTh24 && fUseTrig) title = ";sin^{2}(#theta_{24})";
    else if(paramEnum == cmf::kTh34 && fUseTrig) title = ";sin^{2}(#theta_{34})";
    title += paramScale.second + ";#Delta#chi^{2}";

    TH1F *hist = tfd.make<TH1F>(histName.c_str(),
                                title.c_str(),
                                paramBins.size() - 1,
                                paramBins.data());
    hist->SetMaximum(15.);
    hist->GetXaxis()->CenterTitle();
    hist->GetXaxis()->SetDecimals();

    hist->GetYaxis()->CenterTitle();
    hist->GetYaxis()->SetDecimals();

    // Make TLines to show CLs
    std::vector<Color_t> sigmaColors({kBlue, kRed, kMagenta});
    std::vector<double> cl_percentiles({fContourLevel1Sigma1D,
                                        fContourLevel2Sigma1D,
                                        fContourLevel3Sigma1D});
    std::vector<TLine*> sigmaLines(3);
    for(size_t l = 0; l < sigmaLines.size(); ++l){
     sigmaLines[l]= new TLine(paramBins.front(), cl_percentiles[l], paramBins.back(), cl_percentiles[l]);
     sigmaLines[l]->SetLineStyle(2);
     sigmaLines[l]->SetLineColor(sigmaColors[l]);
    }

    MF_LOG_DEBUG("PlotUtilities")
        << " cl_percentile[0]: "
        << cl_percentiles[0]
        << " "
        << fContourSigmaLabels[0]
        << " cl_percentile[1]: "
        << cl_percentiles[1]
        << " "
        << fContourSigmaLabels[1]
        << " cl_percentile[2]: "
        << cl_percentiles[2]
        << " "
        << fContourSigmaLabels[2];

    // Add a TLegend explain them
    auto *legLines = new TLegend(0.8, 0.1, 0.9, 0.24);
    for(size_t l = 0; l < sigmaLines.size(); ++l){
      legLines->AddEntry(sigmaLines[l], (fContourSigmaLabels[l] + " CL").c_str(), "l");
    }

    legLines->SetBorderSize(0);
    legLines->SetFillStyle(0);

    // make a Canvas to hold it all
    std::string canvName = grName + "Canv";
    auto *canv = tfd.makeAndRegister<TCanvas>(canvName.c_str(),
                                              canvName.c_str(),
                                              2);

    // Log the x axis if this parameter is to be logged - this is a 1D plot
    // so even for the Y parameter we log the x axis
    if     ( xPar && fLogXAxis) canv->SetLogx();
    else if(!xPar && fLogYAxis) canv->SetLogx();

    canv->cd();
    hist->SetStats(false);
    hist->Draw();
    gr->Draw("lsame");
    for(size_t l = 0; l < sigmaLines.size(); ++l) sigmaLines[l]->Draw();
    legLines->Draw();

    // We're done here
    canv->ForceUpdate();
  }

  //......................................................................
  void PlotUtilities::Make2DContours(art::TFileDirectory      & tfd,
                                     std::string         const& chiSqrName,
                                     cmf::PointMap       const& twoDChiSqr,
                                     cmf::GridPoint      const& bestFit,
                                     bool                       useBestFit)
  {
    TH2F *hist = this->Make2DHist(tfd,
                                  chiSqrName,
                                  twoDChiSqr);

    std::vector<double> percentiles({fContourLevel1Sigma2D,
                                     fContourLevel2Sigma2D,
                                     fContourLevel3Sigma2D});
    std::vector<std::string> contourNames({"1sigma",
                                           "2sigma",
                                           "3sigma"});
    std::vector<std::string> contourLabels({fContourSigmaLabels[0],
                                            fContourSigmaLabels[1],
                                            fContourSigmaLabels[2]});


    // check that the minimum for the histogram is less than the lowest contour level,
    // if not, then get rid of that level
    if(hist->GetMinimum() > fContourLevel1Sigma2D){
      percentiles  .erase(percentiles.begin()  );
      contourNames .erase(contourNames.begin() );
      contourLabels.erase(contourLabels.begin());
    }
    

    // MF_LOG_VERBATIM("PlotUtilities")
    //   << "draw contours for ";
    // for(auto const& itr : percentiles)
    //   MF_LOG_VERBATIM("PlotUtilities")
    //  << itr;
    
    hist->SetContour(percentiles.size(), percentiles.data());

    // This avoids disturbing canvases already
    // Yes, we do need a canvas for this to work
    // yes, we do need to call Update() on it
    // No, ROOT is not very friendly.
    TCanvas temp_can;
    temp_can.cd();
    hist->Draw("contlist");
    temp_can.Update();

    auto * plah = (TObjArray *)(gROOT->GetListOfSpecials()->FindObject("contours"));
    if(!plah) return;

    std::vector< std::vector<std::unique_ptr<TGraph> > > graphVec(percentiles.size());

    for(int i = 0; i < plah->GetSize(); ++i){
      std::vector<std::unique_ptr<TGraph> > temp;
      auto * list = (TList*) plah->At(i);
      if(!list) break;
      if(!(list->First())) break;
      for(int igraph = 0; igraph < list->GetSize(); ++igraph){
        if(list->At(igraph)) temp.emplace_back( std::make_unique<TGraph>( *((TGraph*) list->At(igraph))) );
      }
      graphVec[i] = std::move(temp);
    }

    MF_LOG_DEBUG("PlotAssist")
        << "Found : "
        << plah->GetSize()
        << " sets of contours\n"
        << "Parameter 1 and Enum 1: "
        << cmf::cOscParams_Strings[fXParamEnum]
        << " "
        << "Parameter 2 and Enum 2: "
        << cmf::cOscParams_Strings[fYParamEnum]
        << " there are "
        << graphVec.size()
        << " vectors of graphs";

    std::string title = ";";
    title += hist->GetXaxis()->GetTitle();
    title += ";";
    title += hist->GetYaxis()->GetTitle();

    int p = 0;
    int g;
    for(auto const& gritr : graphVec){

      MF_LOG_DEBUG("PlotUtilities")
          << "there are "
          << gritr.size()
          << " graphs for contour "
          << p;

      g = 0;
      for(auto const& itr : gritr){

        std::stringstream name;
        name << std::setprecision(3)
             << chiSqrName
             << "ContourGraph"
             << cmf::cOscParams_Strings[fXParamEnum]
             << cmf::cOscParams_Strings[fYParamEnum]
            << "_CL"
             << contourLabels[p]
             << "_Gr"
             << g;

        MF_LOG_DEBUG("PlotUtilities")
            << "store graph with name "
            << name.str();

        tfd.makeAndRegister<TGraph>(name.str().c_str(),
                                    title.c_str(),
                                    itr->GetN(),
                                    itr->GetX(),
                                    itr->GetY());

        ++g;
      } // end loop over graphs for the percentile
      ++p;
    } // end loop over different percentile graphs

    //MF_LOG_DEBUG << " Found " << graphVec.size() << " standard cl_contour sets\n";

    // Assign colors and line styles
    // have to do some checks because we could have either 2 or 3 sets of graphs
    //-------
    g = 0;
    for(auto & grv : graphVec){
      for(auto & gr : grv){
        if(percentiles.size() == 3 && g == 0){
          gr->SetLineColor(kBlue);
        }
        if((percentiles.size() == 3 && g == 1) || (percentiles.size() == 2 && g == 0)){
          gr->SetLineColor(kRed);
        }
        if((percentiles.size() == 3 && g == 2) || (percentiles.size() == 2 && g == 1)){
          gr->SetLineColor(6);
        }
        gr->SetLineWidth(2);
      }
      ++g;
    }

    // Make the 'backdrop' histogram against which the contours
    // will be drawn

    MF_LOG_DEBUG("PlotUtilities")
        << "Make2DContours function: "
        << "Parameter 1 and Enum 1: "
        << cmf::cOscParams_Strings[fXParamEnum]
        << " "
        << "Parameter 2 and Enum 2: "
        << cmf::cOscParams_Strings[fYParamEnum];

    std::string histName = (chiSqrName                           +
                            cmf::cOscParams_Strings[fXParamEnum] +
                            cmf::cOscParams_Strings[fYParamEnum] +
                            "CLContours");

    TH2F *hBackdrop = tfd.make<TH2F>(histName.c_str(),
                                     title.c_str(),
                                     hist->GetXaxis()->GetNbins(),
                                     hist->GetXaxis()->GetXmin(),
                                     hist->GetXaxis()->GetXmax(),
                                     hist->GetYaxis()->GetNbins(),
                                     hist->GetYaxis()->GetXmin(),
                                     hist->GetYaxis()->GetXmax());

    hBackdrop->GetXaxis()->CenterTitle();
    hBackdrop->GetXaxis()->SetDecimals();

    hBackdrop->GetYaxis()->CenterTitle();
    hBackdrop->GetYaxis()->SetDecimals();
    hBackdrop->GetYaxis()->SetTitleOffset(1.0);

    // Make a new canvas, plot the CL contours,
    //-------
    std::string canName(chiSqrName + "ContoursCanv");
    auto *canv = tfd.make<TCanvas>(canName.c_str(),
                                   (chiSqrName + " Contours").c_str(),
                                   1200,
                                   800);

    hBackdrop->Draw();
    for(auto const& grv : graphVec){
      for(auto const& itr : grv){
        itr->Draw("C SAME");
      }
    }

    // Add a TLegend to describe the CL contours
    // ROOT is inconsistent, and you must add graphs
    // to a TLegend by name, rather than by pointer
    //-------
    TLegend legContours(0.8, 0.1, 0.9, 0.24);
    legContours.SetBorderSize(0);
    legContours.SetFillStyle(0);
    g = 0;
    for(auto & grv : graphVec){
      if(!grv.empty()){
        grv.front()->SetName(contourNames[g].c_str());
        legContours.AddEntry(grv.front().get(),
                             (contourLabels[g] + " CL").c_str(),
                             "l");
      }
      ++g;
    }
    legContours.Draw();

    if(useBestFit){
      // And a TLegend to describe the initial guess and
      // best fit points
      //-------
      TLegend legPoints(0.12, 0.11, 0.27, 0.25);
      legPoints.SetBorderSize(0);
      legPoints.SetFillStyle(0);

      std::stringstream par1_legendentry;
      std::stringstream par2_legendentry;

      std::string Parameter1Label(cmf::cOscParams_Strings_Latex[fXParamEnum]);
      if(fXParamEnum == cmf::kTh23 && fUseTrig) Parameter1Label = "sin^{2}(#theta_{23})";
      if(fXParamEnum == cmf::kTh24 && fUseTrig) Parameter1Label = "sin^{2}(#theta_{24})";

      par1_legendentry
          << Parameter1Label
          << " = "
          << std::setprecision(3)
          << bestFit.X()
          << " "
          << fXParamScale.second;

      std::string Parameter2Label(cmf::cOscParams_Strings_Latex[fYParamEnum]);
      if(fYParamEnum == cmf::kTh23 && fUseTrig) Parameter2Label = "sin^{2}(#theta_{23})";
      if(fYParamEnum == cmf::kTh24 && fUseTrig) Parameter2Label = "sin^{2}(#theta_{24})";

      par2_legendentry
          << Parameter2Label
          << " = "
          << std::setprecision(3)
          << bestFit.Y()
          << " "
          << fYParamScale.second;

      // add the initial and best fit point

      TMarker bestFitMarker(fXParamScale.first * bestFit.X(),
                            fYParamScale.first * bestFit.Y(),
                            8);
      bestFitMarker.SetMarkerColor(1);
      bestFitMarker.Draw();

      legPoints.AddEntry(&bestFitMarker, "Best Fit", "p");
      legPoints.AddEntry(&bestFitMarker, par1_legendentry.str().data(), "");
      legPoints.AddEntry(&bestFitMarker, par2_legendentry.str().data(), "");

      legPoints.Draw();
    }

    if(fLogXAxis) canv->SetLogx();
    if(fLogYAxis) canv->SetLogy();

    canv->Update();
    canv->Write();
  }

  //......................................................................
  void PlotUtilities::MakeCLHeatMaps(art::TFileDirectory               & tfd,
                                     std::string                  const& baseName,
                                     std::vector<cmf::ChiSqrInfo> const& chiSqrInfoVec)
  {

    // make a 2D histogram for each CL
    // make a 2D histogram and fill it with the Delta chiSqr values
    std::string name = (baseName                             +
                        cmf::cOscParams_Strings[fXParamEnum] +
                        cmf::cOscParams_Strings[fYParamEnum] +
                        "1sigmaHeatMap");

    std::vector<TH2F*> hmVec(3);

    for(size_t i = 0; i < hmVec.size(); ++i){

      if(i > 0)
        name.replace(name.find(std::to_string(i) + "sigma"), 1, std::to_string(i+1));

      hmVec[i] = this->MakeGeneric2DHist(tfd,
                                         name);
      hmVec[i]->SetTitle((std::to_string(i + 1) + "#sigma").c_str());
    }

    // loop over the chiSqrInfo structs to find the delta chiSqr at each point
    MF_LOG_DEBUG("PlotUtilities")
      << "there are "
      << chiSqrInfoVec.size()
      << " chiSqrInfo objects to loop over";

    size_t uni = 0;
    for(auto const& csItr : chiSqrInfoVec){
      for(auto const& twoDItr : csItr.twoDChiSqr){

        MF_LOG_DEBUG("PlotUtilities")
          << "twoDChiSqr for universe "
          << uni
          << " "
          << twoDItr.first
          << " "
          << twoDItr.second;
        
        if(1.0-TMath::Prob(twoDItr.second, 2) < fContourLevel1Sigma2D)
          hmVec[0]->Fill(twoDItr.first.X() * fXParamScale.first,
                         twoDItr.first.Y() * fYParamScale.first);
        if(1.0-TMath::Prob(twoDItr.second, 2) < fContourLevel2Sigma2D)
          hmVec[1]->Fill(twoDItr.first.X() * fXParamScale.first,
                         twoDItr.first.Y() * fYParamScale.first);
        if(1.0-TMath::Prob(twoDItr.second, 2) < fContourLevel3Sigma2D)
          hmVec[2]->Fill(twoDItr.first.X() * fXParamScale.first,
                         twoDItr.first.Y() * fYParamScale.first);
      }
      ++uni;
    } // end loop over ChiSqrInfos

    // make a canvas to hold these histograms
    name = baseName                             +
           cmf::cOscParams_Strings[fXParamEnum] +
           cmf::cOscParams_Strings[fYParamEnum] +
           "HeatMapCanv";

    auto *canv = tfd.make<TCanvas>(name.c_str(),
                                   name.c_str(),
                                   800,
                                   800);

    if(fLogXAxis) canv->SetLogx();
    if(fLogYAxis) canv->SetLogy();
    canv->Divide(2,2);

    for(size_t i = 0; i < hmVec.size(); ++i){
      canv->cd(i + 1);
      hmVec[i]->Draw("colz");
    }

    canv->Update();
    canv->Write();

    return;
  }

  //......................................................................
  // Normalize the bin contents to be the equivalent to the desired width
  void PlotUtilities::NormalizeBinContents(TH1*   hist,
                                           double normToVal)
  {
    if(!fNormalizeBins) return;

    double norm;
    for(int b = 0; b < hist->GetNbinsX(); ++b){
      norm = normToVal / hist->GetBinWidth(b + 1);
      hist->SetBinContent(b + 1, hist->GetBinContent(b + 1) * norm);
      hist->SetBinError(b + 1, hist->GetBinError(b + 1) * norm);
    } // end loop over bins
  }

  //......................................................................
  void PlotUtilities::MakeEnergySpectraFromBins(cmf::Spectrum         const& binSpec,
                                                cmf::Spectrum         const& binCount,
                                                std::map<long, TH1D*>      & energySpecMap,
                                                std::string           const& name,
                                                art::TFileDirectory        & tfd,
                                                std::string           const& uniqueID)
  {
    std::string   histName;
    std::string   histTitle;
    cmf::Spectrum bins;
    TH1D *hist;

    // loop over the keys we are using to determine
    // the offset for each selection
    long offset;

    double binE;
    double scale;

    cmf::SelectionType_t sel;
    cmf::DetType_t       det;
    cmf::BeamType_t      beam;

    for(auto const& itr : cmf::SelectionUtility::Instance()->SelectionsToUseAsKeys()){
      offset = cmf::CovarianceBinUtility::Instance()->KeyToOffset(itr);
      sel    = cmf::KeyToSelectionType(itr);
      det    = cmf::KeyToDetectorType(itr);
      beam   = cmf::KeyToBeamType(itr);
      scale = (det == cmf::kNEARDET) ? 1.e-3 : 1.;

      bins.clear();
      cmf::CovarianceBinUtility::Instance()->SelectionHistBinning(sel,
                                                                  det,
                                                                  bins);

      // check that we have a sane number of bins
      if(offset + cmf::CovarianceBinUtility::Instance()->SelectionHighEdges(itr).size() > binSpec.size())
        throw cet::exception("CovarianceFitHelper")
            << "Filling data energy spectra appears to want more bins"
            << " than are in the data bin spectrum: "
            << offset + cmf::CovarianceBinUtility::Instance()->SelectionHighEdges(itr).size()
            << " vs "
            << binSpec.size();

      histName = (name +
                  uniqueID +
                  cmf::cBeamType_Strings[beam] +
                  cmf::cDetType_Strings[det]   +
                  cmf::cSelectionType_Strings[sel]);

      histTitle = (cmf::cBeamType_LatexStrings[beam]     + " " +
                   cmf::cDetType_Strings[det]            + " " +
                   cmf::cSelectionType_LatexStrings[sel] + ";Energy (GeV);");

      histTitle += (scale < 1) ? "#times10^{3} Events" : "Events";

      if(name.find("MC") != std::string::npos) histTitle.insert(0, "MC ");

      MF_LOG_DEBUG("CovarianceFitHelper")
          << "total Energy spectrum hist is "
          << histName
          << " "
          << itr;

//    for(auto const& binitr : bins)
//      LOG_VERBATIM("CovarianceFitHelper")
//          << "bin edge: "
//          << binitr;

      if(energySpecMap.count(itr) < 1)
        hist = energySpecMap.emplace(itr, tfd.make<TH1D>(histName.c_str(),
                                                         histTitle.c_str(),
                                                         bins.size() - 1,
                                                         bins.data())).first->second;
      else
        hist = energySpecMap.find(itr)->second;

      hist->GetXaxis()->CenterTitle();
      hist->GetYaxis()->CenterTitle();
      if(name.find("MC") != std::string::npos){
        hist->SetMarkerColor(kRed);
        hist->SetLineColor(kRed);
      }

      for(size_t b = offset; b < offset + cmf::CovarianceBinUtility::Instance()->SelectionHighEdges(itr).size(); ++b){
        binE = cmf::CovarianceBinUtility::Instance()->BinToEnergy(b);
        hist->Fill(binE, binSpec[b] * scale);

        // the uncertainty for each bin is the fractional uncertainty for the total count in the
        // bin (1 / sqrt(binCount[b]) multiplied by the total weight in the bin (binSpec[b])
        // only do it if the bin count > 0
        if(binCount[b] > 0.)
          hist->SetBinError(hist->FindBin(binE),
                            scale * (binSpec[b] / std::sqrt(binCount[b])));
      }// end filling the spectrum

      MF_LOG_DEBUG("PlotUtilities") 
        << "histogram " 
        << hist->GetName() 
        << " has " 
        << hist->Integral() 
        << " selected events";

      // now normalize the bin contents if desired to the desired width
      if(cmf::IsNuESelected(sel))
        this->NormalizeBinContents(hist, fNuENormVal);
      else if(cmf::IsNuMuSelected(sel))
        this->NormalizeBinContents(hist, fNuMuNormVal);
      else if (cmf::IsNCSelected(sel))
        this->NormalizeBinContents(hist, fNCNormVal);

    } // end loop over selections being used in the analysis
  }

  //......................................................................
  TH1F* PlotUtilities::MakeSpectrumHistogram(art::TFileDirectory      & tfd,
                                             std::string         const& baseName,
                                             std::string         const& baseTitle,
                                             cmf::Spectrum       const& spectrum)
  {
    if(spectrum.size() < 1){
      MF_LOG_VERBATIM("PlotUtilities")
          << "spectrum size for "
          << baseName
          << " is 0, are you sure you filled it?";
      return nullptr;
    }
    
    auto spectrumHist = tfd.make<TH1F>(baseName.c_str(),
                                       baseTitle.c_str(),
                                       spectrum.size(),
                                       0, 
                                       spectrum.size());

    spectrumHist->Sumw2();

    for(size_t b = 0; b < spectrum.size(); ++b) spectrumHist->Fill(b, spectrum[b]);

    return spectrumHist;
  }

  //......................................................................
  // The returned histogram has the axes and axis ranges defined by fXParam(2)
  // only for now
  TH2F* PlotUtilities::MakeGeneric2DHist(art::TFileDirectory      & tfd,
                                         std::string         const& histName)
  {
    // Delta m^2 is always on the y axis for numu
    // angle is always on y axis for nue contours

    std::string XaxisTitle = cmf::cOscParams_Strings_Latex[fXParamEnum] + fXParamScale.second;
    std::string YaxisTitle = cmf::cOscParams_Strings_Latex[fYParamEnum] + fYParamScale.second;

    if(cmf::IsAngleParameter(fXParamEnum) && fUseTrig) XaxisTitle = "sin^{2}(" + cmf::cOscParams_Strings_Latex[fXParamEnum] + ")";
    if(cmf::IsAngleParameter(fYParamEnum) && fUseTrig) YaxisTitle = "sin^{2}(" + cmf::cOscParams_Strings_Latex[fYParamEnum] + ")";

    std::string title = (";" + XaxisTitle +
                         ";" + YaxisTitle);


    MF_LOG_DEBUG("ContourMaker")
        << "Creating 2D histogram with range "
        << fXParamBinEdges.front() << " " << fXParamBinEdges.back() << " "
        << fYParamBinEdges.front() << " " << fYParamBinEdges.back();

    auto* hist2D =  tfd.make<TH2F>(histName.c_str(),
                                   title.c_str(),
                                   fXParamDivs+1,
                                   fXParamBinEdges.data(),
                                   fYParamDivs+1,
                                   fYParamBinEdges.data());

    hist2D->GetXaxis()->CenterTitle();
    hist2D->GetXaxis()->SetDecimals();
    hist2D->GetYaxis()->CenterTitle();
    hist2D->GetYaxis()->SetDecimals();

    return hist2D;
  }

  //......................................................................
  TH2F* PlotUtilities::Make2DHist(art::TFileDirectory      & tfd,
                                  std::string         const& chiSqrName,
                                  cmf::PointMap       const& twoDChiSqr)
  {
    // make a 2D histogram and fill it with the Delta chiSqr values
    std::string histName = (chiSqrName                           +
                            cmf::cOscParams_Strings[fXParamEnum] +
                            cmf::cOscParams_Strings[fYParamEnum] +
                            "DeltaChiSqr");

    auto* hist_deltachisq = this->MakeGeneric2DHist(tfd,
                                                    histName);

    histName = (chiSqrName                           +
                cmf::cOscParams_Strings[fXParamEnum] +
                cmf::cOscParams_Strings[fYParamEnum] +
                "DeltaChiSqr_Prob");

    auto* hist = this->MakeGeneric2DHist(tfd,
                                         histName);

    double dcsq;
    for(auto const& itr : twoDChiSqr){

      MF_LOG_DEBUG("ContourMaker")
          << itr.first.X()
          << " "
          << itr.first.Y()
          << " "
          << itr.second;

      dcsq = (itr.second < 0.1) ? 0.1 : itr.second;
      hist_deltachisq->Fill(itr.first.X() * fXParamScale.first,
                            itr.first.Y() * fYParamScale.first,
                            dcsq);

      hist->Fill(itr.first.X() * fXParamScale.first,
                 itr.first.Y() * fYParamScale.first,
                 1.0-TMath::Prob(itr.second, 2));
    }

//    for(int b = 0; b < hist_deltachisq->GetXaxis()->GetNbins(); ++ b)
//      MF_LOG_VERBATIM("CMFContourMaker")
//      << "x axis "
//      << b + 1
//      << " "
//      << hist_deltachisq->GetXaxis()->GetBinCenter(b + 1);
//
//    for(int b = 0; b < hist_deltachisq->GetYaxis()->GetNbins(); ++ b)
//      MF_LOG_VERBATIM("CMFContourMaker")
//          << "y axis "
//          << b + 1
//          << " "
//          << hist_deltachisq->GetYaxis()->GetBinCenter(b + 1);

    return hist;
  }

  //......................................................................
  void PlotUtilities::Make2DHiddenParHistogram(art::TFileDirectory      & tfd,
                                               std::string         const& baseName,
                                               cmf::PointMap       const& parVals)
  {
    if(parVals.size() < 1){
      MF_LOG_VERBATIM("PlotUtilities")
          << "point map size for "
          << baseName
          << " is 0, are you sure you filled it?";
      return;
    }

    auto* parValHist = this->MakeGeneric2DHist(tfd,
                                               baseName);

    for(auto const& itr : parVals){
      parValHist->Fill(itr.first.X() * fXParamScale.first,
                       itr.first.Y() * fYParamScale.first,
                       itr.second);
    }
  }

} // end cmf namespace