KinematicsPlots.C

Go to the documentation of this file.

#ifdef __CINT__
void KinematicsPlots(std::string hornCurrent, std::string prodName)
{
  std::cout << "CINT not supported (or tolerated) here" << std::endl;
}
#else

#include <iostream>
#include <map>
#include <memory>
#include <string>
#include <typeinfo>
#include <unordered_map>
#include <unordered_set>

#include "boost/algorithm/string.hpp"
#include "boost/tokenizer.hpp"

#include "TCanvas.h"
#include "TError.h"
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "THStack.h"
#include "TKey.h"
#include "TLegend.h"
#include "TLegendEntry.h"
#include "TList.h"
#include "TString.h"
#include "TLatex.h"
#include "TROOT.h"

#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/Utilities.h"

// /////////////////////////////////////
//
//  custom types for use below

const std::map<std::string, int> INT_COLORS = {
        { "QE",    kAzure+7 },
        { "RES",   kGreen+2 },
        { "DIS",   kGray },
        { "MEC",   kOrange-1 },
        { "Other", kBlack }
};

const std::vector<std::string> INT_ORDER = {
        "Other",
        "DIS",
        "RES",
        "QE",
        "MEC"
};

const std::map<std::string, std::string> INT_LEGEND = {
        { "QE",    "QE" },
        { "RES",   "RES" },
        { "DIS",   "DIS" },
        { "MEC",   "NOvA 2p2h" },
        { "Other", "Other" }
};

// would probably be better to determine all of these dynamically from the keys in the file...
// but I don't feel like it right now
const std::vector<std::string> MECQE_WGTS = {
    "CV2018",
    "CV2018-NoMEC",
    "NOvA-MECdown",
    "NOvA-MECup",
    "MINERvA-1p1h",
    "MINERvA-2p2h",
    "MINERvA-2p2hnp",
    "MINERvA-2p2hpp",
};

// again, would be nice to determine dynamically, but this is easier
const std::unordered_map<std::string, double> Q0Q3_MAX = {
    {"FHC", 8e4},
//    {"FHC", 2e3},
    {"RHC", 3e4},
};

// scale down by this much.  hard-code for now...
const double Q0Q3_SCALE = 1e-4;

const double MC_NOMINAL_POT = 9e20;  // for MC-only plots

const double Q3_SLICE_SIZE = 0.1;
const unsigned int Q0Q3_SLICE_OFFSET = 1;

class PlotKey
{
        public:
                PlotKey() = default;
                PlotKey(const std::string & wgtnm, const std::string & intnm)
                        : _mecWgt(wgtnm), _interactionType(intnm),
                          _combination(Form("{%s}{%s}", wgtnm.c_str(), intnm.c_str()))
    {};

                PlotKey(const PlotKey & other) = default;
                        
                bool operator<(const PlotKey& other) const {
                  if (this->_mecWgt == other.mecWgt())
                    return this->_interactionType < other.interactionType();
                  else
                    return this->_mecWgt < other.mecWgt();
                };

                std::string mecWgt() const { return this->_mecWgt; };
                std::string interactionType() const { return this->_interactionType; };
                
                const std::string & combination() const { return this->_combination; };
        private:
                std::string _mecWgt;
                std::string _interactionType;
        
                // need this for ordering in the map... sigh.
                std::string _combination;
};

std::ostream & operator<<(std::ostream &os, PlotKey const &pk)
{
  return os << pk.combination();
}

// /////////////////////////////////////
//
//  helper obj


class Plotter
{
  public:
    Plotter(const std::string & hornCurrent, const std::string & prodName, const std::string & outDir)
      : fHornCurrent(hornCurrent), fProdName(prodName),
        fOutDir( Form("%s/%s/%s", outDir.c_str(), hornCurrent.c_str(), prodName.c_str()) ),
        fHistFileData(TFile::Open( Form("%s/%s/rw_histos_%s_data.root", outDir.c_str(), hornCurrent.c_str(), fProdName.c_str()), "recreate" )),
        fHistFileMC(TFile::Open( Form("%s/%s/rw_histos_%s_MC.root", outDir.c_str(), hornCurrent.c_str(), fProdName.c_str()), "recreate" ))
    {
      this->GetDataPlots();
      this->GetMCPlots();
    }

    ~Plotter()
    {
      fHistFileData->Close();
      fHistFileMC->Close();

      delete fHistFileData;
      delete fHistFileMC;
    }

  private:
    std::string fHornCurrent;
    std::string fProdName;
    std::string fOutDir;

    std::map<std::string, std::unique_ptr<ana::Spectrum>> fDataPlots;
    std::map<std::string, std::map<PlotKey,std::unique_ptr<ana::Spectrum>>> fMCPlots;

    std::map<PlotKey, std::vector<std::unique_ptr<TCanvas>>> fq0q3Canvases;  //ugh

    TFile * fHistFileData;
    TFile * fHistFileMC;

    std::unordered_set<std::string> fdataPlotsWritten;
    std::unordered_set<std::string> fMCPlotsWritten;

    void WritePlotToHistFile(const std::string & plotName, TH1 * hist, bool isData=true)
    {
      auto name = plotName.c_str();

      TFile * f = (isData) ? fHistFileData : fHistFileMC;
      decltype(fdataPlotsWritten) & cache = isData ? fdataPlotsWritten : fMCPlotsWritten;

      if (cache.find(name) != cache.end())
        return;

      f->cd();
      hist->SetName(name);
      hist->Write(name);
      cache.insert(name);
    }

    void GetDataPlots()
    {
      ana::DontAddDirectory inclGuard;
      TFile f( Form("/nova/ana/users/jwolcott/2p2h_tuning/2018/hists/rwgt_spectra_%s_%s_Data.root", fProdName.c_str(), fHornCurrent.c_str()) );
      TList * keys = f.GetListOfKeys();

      fDataPlots.clear();
      TIter it(keys);
      std::for_each(it.Begin(), TIter::End(), [this](TObject * listItem){
        TKey * key = dynamic_cast<TKey*>(listItem);

        TDirectory * dir = dynamic_cast<TDirectory*>(key->ReadObj());
        if (!dir)
          return;

        fDataPlots.emplace(key->GetName(), std::move(ana::Spectrum::LoadFrom(dir)));
      });
    }

    ////////

    bool SplitMCPlotName(const std::string& keyedName, std::string & name, PlotKey & plotkey)
    {
      std::vector<std::string> fields;

      typedef boost::tokenizer<boost::char_separator<char> >
          tokenizer;
      boost::char_separator<char> sep("{}");
      tokenizer tokens(keyedName, sep);

//      std::cout << "Parsing keyed name: " << keyedName << std::endl;
      for (const auto & f : tokens)
      {
    //    std::cout << "  raw = " << f;
        if (std::count(f.begin(), f.end(), '=') < 2)  // if it's a named category, it has a 'group=' and a 'cat='
          fields.push_back(f);
        else
        {
          std::vector<std::string> subfields;
          boost::split(subfields, f, boost::is_any_of("="), boost::token_compress_on);
          fields.push_back(subfields.back());
        }

    //    std::cout << ";  parsed = " << fields.back() << std::endl;
      }

      // truth variables have only a var name and an interaction class...
      if (fields.size() == 2)
        plotkey = PlotKey("", fields[1]);
      // ... while other plots have all the relevant fields: mec weight name, interaction class...
      else if (fields.size() == 3)
        plotkey = PlotKey(fields[1], fields[2]);
      else
      {
        std::cerr << Form("Warning: Can't parse Spectrum name: '%s'", keyedName.c_str()) << std::endl;
        return false;
      }

      name = fields[0];

//       std::cout << "  name, computed key = " << name << ", " << plotkey << std::endl;

      return true;
    }

    ////////

    void GetMCPlots()
    {
      ana::DontAddDirectory inclGuard;
      fMCPlots.clear();
      auto fileName = Form("/nova/ana/users/jwolcott/2p2h_tuning/2018/hists/rwgt_spectra_%s_%s_MC.root",
                           this->fProdName.c_str(), this->fHornCurrent.c_str());
      TFile f(fileName);
      if (f.IsZombie())
        return;

      TList * keys = f.GetListOfKeys();

      TIter it(keys);
      std::for_each(it.Begin(), TIter::End(), [this](TObject * listItem){
        TKey * key = dynamic_cast<TKey*>(listItem);

        TDirectory * dir = dynamic_cast<TDirectory*>(key->ReadObj());
        if (!dir)
          return;

        std::string name;
        PlotKey plotkey;
        bool success = SplitMCPlotName(key->GetName(), name, plotkey);
        if (!success)
          return;

        // if slot not already occupied OR preceding entry has nothing in it
        if (fMCPlots[name].find(plotkey) == fMCPlots[name].end() || fMCPlots[name].at(plotkey)->ToTH1(fMCPlots[name].at(plotkey)->POT())->GetEntries() == 0)
          fMCPlots[name][plotkey] = std::move(ana::Spectrum::LoadFrom(dir));
      });

      std::cout << "MC plots: ";
      for (const auto & plotPair : fMCPlots)
      {
        std::cout << plotPair.first << ":\n";
    //    for (const auto & keyPair : plotPair.second)
    //      std::cout << "   " << keyPair.first << "\n";
        std::cout << std::endl;
      }
    //  throw std::runtime_error("look these over first");

    }

    ////////

    TCanvas * DataMCComparison(const std::string & plotName,
                               const PlotKey & key,
                               ana::Spectrum * dataPlot,
                               const std::map<PlotKey,
                               const ana::Spectrum*> & mcplots,
                               const PlotKey & plotKey
    )
    {
      // these guys are unfortunately going to get leaked.
      // it's too much work otherwise to make sure they persist
      // through this function's return until the Print()ing
      // of the canvas.  this macro is small anyway
      THStack * stack = new THStack;
      stack->SetName("stack");
      TLegend * legend = new TLegend(0.7, 0.5, 0.9, 0.85);
      legend->SetFillStyle(0);

      // this one, on the other hand, is returned and stuck
      // into a unique_ptr, so he's ok
      TCanvas * canvas = new TCanvas;
      double dataPOT = dataPlot->POT();
      static bool announcedPOT = false;
      if (!announcedPOT)
      {
         std::cout << " data POT = " << dataPOT << std::endl;
         announcedPOT = true;
      }
      TH1 * dataTH1 = dataPlot->ToTH1(dataPOT);
      dataTH1->SetLineWidth(2);
      dataTH1->SetMarkerStyle(8);
      if (std::string(dataTH1->GetXaxis()->GetTitle()).find("\"q0\"") != std::string::npos)  // this is a ridiculous way to decide if this is one of the q0,q3 slices, but the histogram name isn't available here
        dataTH1->SetMarkerSize(1.2);
      else
        dataTH1->SetMarkerSize(0.8);
      dataTH1->Draw("pe");  // draw the data first so that the labeled axes show up.  don't feel like copying them to the THStack
      this->WritePlotToHistFile(plotName.c_str(), dataTH1, true);

      for (const auto & intClass : INT_ORDER)
      {
        PlotKey key(plotKey.mecWgt(), intClass);
        TH1 * mcTH1 = mcplots.at(key)->ToTH1(dataPOT);
        mcTH1->SetLineWidth(1);
//        std::cout << " integral = " << mcTH1->Integral() << std::endl;
        mcTH1->SetFillColor(INT_COLORS.at(intClass));
        this->WritePlotToHistFile(plotName + key.combination(), mcTH1, false);
        stack->Add(mcTH1, "hist");
        // do this the hard way to get the legend order the same as the stack instead of reversed
        if (plotKey.mecWgt().find("NoMEC") == std::string::npos || intClass != "MEC")
           legend->GetListOfPrimitives()->AddFirst(new TLegendEntry(mcTH1, INT_LEGEND.at(intClass).c_str(), "f"));
      }
      stack->Draw("same");

      TH1 * mcsum = dynamic_cast<TH1*>((*stack->GetStack())[stack->GetStack()->LastIndex()]);

      auto maxVal = std::max(dataTH1->GetMaximum(), mcsum->GetMaximum());

      maxVal *= 1.2;
      dataTH1->SetMaximum(maxVal);
      dataTH1->GetYaxis()->SetRange(0, maxVal);

      dataTH1->Draw("pe same"); // yep, do it again, to move the points back up to the top

      legend->Draw();

      return canvas;
    }

    ////////

  public:

    void DataMCComparisons()
    {
      for (const auto & dataPair : fDataPlots)
      {
        const std::string & plotName = dataPair.first;
        const std::unique_ptr<ana::Spectrum> & dataPlot = dataPair.second;

        for (const auto & mecWgtName : MECQE_WGTS)
        {
            if (fMCPlots.find(plotName) == fMCPlots.end())
            {
              std::cerr << Form("  warning: MC plot with mec wgt '%s' not found for data plot '%s'", mecWgtName.c_str(), plotName.c_str()) << std::endl;
              continue;
            }

    //                          std::cout << "  MC plots " << (MCPlots.find(plotName) == MCPlots.end() ? "DO NOT have" : "HAVE") << " plot with key '" << plotName << "' in them" << std::endl;
            PlotKey key(mecWgtName, "");
            const auto & plotCollection = fMCPlots.at(plotName);

    //                          std::cout << "  MC plots with name '" << plotName << "' are: " << std::endl;
    //                          for (const auto & plotPair : plotCollection)
    //                            std::cout << "     " << plotPair.first << std::endl;;

              // find the set of keys with the same weights but possibly different interaction types
    //                            typename std::remove_reference<decltype(cont)>::type
              std::map<PlotKey, const ana::Spectrum*> plotsByInteraction;
              for (const auto & plotPair : plotCollection)
              {
    //                              std::cout << "      candidate: " << plotPair.first;
                if ( plotPair.first.mecWgt() != key.mecWgt())
                {
    //                                std::cout << " ... NOPE" << std::endl;
                  continue;
                }
    //                              std::cout << "  YES" << std::endl;
                plotsByInteraction.insert( std::make_pair(plotPair.first, plotPair.second.get()) );
//                std::cout << plotPair.first << " " << plotPair.second->Integral(plotPair.second->POT()) << std::endl;;
              }

              try
              {
                std::unique_ptr<TCanvas> dataMCComparison(DataMCComparison(plotName, key, dataPlot.get(), plotsByInteraction, key));

                std::size_t q3Loc = plotName.find("q3=");
                if (q3Loc != std::string::npos)
                {
                  TIter iter(dataMCComparison->GetListOfPrimitives());
                  TObject * obj;
                  while ( (obj = iter()) )
                  {
                    TH1 * h = dynamic_cast<TH1*>(obj);
                    if (!h)
                      continue;
                    h->GetXaxis()->SetRangeUser(0, 1);
                    h->GetYaxis()->SetRangeUser(0, Q0Q3_MAX.at(this->fHornCurrent));
                  }
                }
                for (const auto & ext : {"png", "root", "eps", "pdf"})
                  dataMCComparison->Print( Form("%s/%s/%s.%s", fOutDir.c_str(), mecWgtName.c_str(), plotName.c_str(), ext) );

                // stash the q0q3 plots away for later
                if (q3Loc != std::string::npos)
                {
                  auto lb = std::atof(plotName.substr(q3Loc+3, plotName.find("-")).c_str());
                  auto & canvasColl = this->fq0q3Canvases[key];
                  int q3SliceIdx = 10*lb;
                  if (q3SliceIdx > (int)(canvasColl.size())-1)  // gotta watch out for tricky integer overflows
                    canvasColl.resize(q3SliceIdx+1);
//                  std::cout << " for plot " << plotName << ", concluded q3 start = " << lb << ", slice number = " << q3SliceIdx << std::endl;
                  canvasColl[q3SliceIdx] = std::move(dataMCComparison);
                }
              }
              catch (std::out_of_range & e)
              {
                std::cout << "  couldn't find plot '" << plotName << "' with key '" << key << "'.  Skipping data-MC comparison." << std::endl;
                continue;
              }
        } // for (mecWgtName)
      } // for (dataPair)
    } // DataMCComparisons()

    ////////

    void q0q3Panel()
    {
      if (this->fq0q3Canvases.empty())
        return;

      double q0max = Q0Q3_MAX.at(this->fHornCurrent) * Q0Q3_SCALE;
      q0max *= 0.98; // for axis autolabeling purposes

      for (const auto & keyPair : this->fq0q3Canvases)
      {
        auto & key = keyPair.first;
        auto & canvases = keyPair.second;
        std::cout << "  making q0q3 panel for key: " << key << std::endl;

        const double bottomMarginMaster = 0.12;
        const double leftMarginMaster = 0.1;
        const double rightMarginMaster = 0.05;
        const double topMarginMaster = 0.05;

        TCanvas masterC("q0q3panel", "(q0, |q|) panel", 700, 700);
        auto nCanvases = canvases.size();
//        masterC.Divide(3, nCanvases / 3, 0, 0);
                auto nRows = nCanvases / 3;
                double padWidth = (1. - leftMarginMaster - rightMarginMaster) / 3;
                double padHeight = (1. - topMarginMaster - bottomMarginMaster) / nRows;

        // yes, skip 0.  nothing interesting in that bin anywho
        for (unsigned int cIdx = Q0Q3_SLICE_OFFSET; cIdx < nCanvases; cIdx++)
        {
          if (!canvases[cIdx])
          {
            std::cout << "   null" << std::endl;
            continue;
          }

          TIter it_prim(canvases[cIdx]->GetListOfPrimitives());
          bool adjustedData = false;  // the data histogram sometimes is in there twice
          while (TObject * obj = it_prim())
          {
            if (auto l = dynamic_cast<TLegend*>(obj))
                  canvases[cIdx]->GetListOfPrimitives()->Remove(l);
            if ( auto h = dynamic_cast<TH1*>(obj) )
            {

              h->SetTitle(";;");
              if (!adjustedData)
              {
                      adjustedData = true;
                      h->Scale(Q0Q3_SCALE);
              }
              h->GetYaxis()->SetRangeUser(0, q0max);
              h->GetYaxis()->SetDecimals();

              if ((cIdx - Q0Q3_SLICE_OFFSET) % 3 == 0)
              {
                h->SetLabelSize(0.03, "Y");
                h->GetYaxis()->ChangeLabel(-1, 0, 0);  // erase the last axis label so it doesn't clash with the 0 above it
              }
              else
                 h->SetLabelSize(0, "Y");

              if ((cIdx-Q0Q3_SLICE_OFFSET) / 3 >= (nCanvases-Q0Q3_SLICE_OFFSET) / 3 - 1)
              {
                    h->SetLabelSize(0.03, "X");
                    h->GetXaxis()->ChangeLabel(-1, -1, 0);  // erase the last axis label so it doesn't clash with the 0 next to it
              }
              else
                h->SetLabelSize(0, "X");

              h->GetXaxis()->SetRangeUser(0, 0.48); // not 0.5 so the last label isn't shown
              h->GetXaxis()->SetDecimals();
            }
            else if ( auto st = dynamic_cast<THStack*>(obj) )
                {
                      for (auto o : *st->GetHists())
                      {
                        auto h = dynamic_cast<TH1*>(o);
                            h->Scale(Q0Q3_SCALE);
                            h->GetYaxis()->SetRangeUser(0, q0max);
                      }
                      st->Modified();
                }
          }

          auto idx = cIdx - Q0Q3_SLICE_OFFSET;
          auto column = idx % 3;
          double leftMargin = leftMarginMaster + (column * padWidth);
          double rightMargin = rightMarginMaster + (3 - column - 1) * padWidth;
          auto row = idx / 3;
          double topMargin = topMarginMaster + row * padHeight;
          double bottomMargin = bottomMarginMaster + (nRows - row - 1) * padHeight;

          masterC.cd();
          TLatex * l = new TLatex(0.25, q0max*0.9, Form("%.1f #leq |#vec{q}|/(GeV/c) < %.1f", cIdx * Q3_SLICE_SIZE, (cIdx+1)*Q3_SLICE_SIZE));
          l->SetTextAlign(22); // middle-middle
          l->SetTextSize(0.025);
          canvases[cIdx]->GetListOfPrimitives()->Add(l);
//          canvases[cIdx]->GetListOfPrimitives()->Print();

          canvases[cIdx]->Update();

          auto p = new TPad(std::to_string(cIdx).c_str(), "", 0, 0, 1, 1);
          p->cd();
//          gROOT->SetSelectedPad(canvases[cIdx].get());
          auto c = canvases[cIdx]->DrawClonePad();
//          gROOT->GetListOfCanvases()->Remove(c);
                p->SetMargin(leftMargin, rightMargin, bottomMargin, topMargin);
//              std::cout << " left right top bottom margins are " << leftMargin << " " << rightMargin << " " << topMargin << " " << bottomMargin << std::endl;
                p->SetFillStyle(0);
                p->RedrawAxis();

                masterC.cd();
          p->Draw();

 //         canvases[cIdx]->Print();

        }

            masterC.cd();

        // rebuild the legend now
        TLegend leg(leftMarginMaster + 0.3*padWidth,  1 - (topMarginMaster+0.95*padHeight),
                    leftMarginMaster + padWidth,      1 - (topMarginMaster+0.2*padHeight));
        leg.SetFillStyle(0);
        auto pad = dynamic_cast<TPad*>(masterC.GetPrimitive("1"));
        auto stack = dynamic_cast<THStack*>(pad->GetPrimitive("stack"));
        for (auto o : *stack->GetHists())
        {
                auto h = dynamic_cast<TH1*>(o);
                auto toks = TString(h->GetName()).Tokenize("{}");
                leg.GetListOfPrimitives()->AddFirst(new TLegendEntry(h, INT_LEGEND.at(toks->Last()->GetName()).c_str(), "f"));
        }
        TH1 * h;
        for (auto o : *pad->GetListOfPrimitives())
        {
                if ( (h = dynamic_cast<TH1*>(o)) )
                        break;
        }
        leg.GetListOfPrimitives()->AddFirst(new TLegendEntry(h, "ND Data", "lpe"));
        leg.Draw();

        TLatex xaxisTitle(0.5, 0.03, "Visible E_{had} (GeV)");
        xaxisTitle.SetNDC();
        xaxisTitle.SetTextAlign(22);
        xaxisTitle.SetTextSize(0.04);
        xaxisTitle.Draw();
        TLatex yaxisTitle(0.025, 0.5, Form("10^{%d} Events", int(-log10(Q0Q3_SCALE))));
        yaxisTitle.SetTextAngle(90);
        yaxisTitle.SetNDC();
        yaxisTitle.SetTextAlign(22);
        yaxisTitle.SetTextSize(0.04);
        yaxisTitle.Draw();

        TLatex beamText(leftMarginMaster + 0.02, 1 - topMarginMaster + 0.01,
                        (fHornCurrent == "FHC" ? "Neutrino beam" : "Antineutrino beam"));
        beamText.SetNDC();
        beamText.SetTextAlign(11);  // align left-bottom
        beamText.SetTextSize(0.03);
        beamText.Draw();

        for (const auto & ext : {"png", "root", "eps", "pdf"})
          masterC.Print(Form("%s/%s/%s.%s", fOutDir.c_str(), key.mecWgt().c_str(), "q0q3panel", ext));
      }
    }

    ////////

    void TruthPlots()
    {
      TCanvas c;
      // migration matrices: don't stack these
        for (const auto & intName : INT_ORDER)
        {
          c.Clear();
          if (fMCPlots.find("Q3Migration") != fMCPlots.end())
          {
          auto & spectrum = fMCPlots.at("Q3Migration").at(PlotKey("NoMECWgt", intName));
          TH2 * mcTH2 = spectrum->ToTH2(MC_NOMINAL_POT);
          mcTH2->Draw("colz");
          for (const auto & ext : {"png", "root", "eps", "pdf"})
            c.Print( Form("%s/NoMECWgt/Q3Migration_%s.%s", fOutDir.c_str(), intName.c_str(), ext) );
          }

          // these only exist for QE & Dytman-MEC
          if (intName != "QE" && intName != "MEC")
            continue;
          if (fMCPlots.find("TrueQ0VsTrueQmag") == fMCPlots.end())
                continue;
          for (const auto & mecWgtName : MECQE_WGTS)
          {
            c.Clear();
            PlotKey pk(mecWgtName, intName);
//            std::cout << "  trying plot: TrueQ0VsTrueQmag/" << pk << std::endl;
            auto & spectrum = fMCPlots.at("TrueQ0VsTrueQmag").at(PlotKey(pk));
            TH2 * mcTH2 = spectrum->ToTH2(MC_NOMINAL_POT);
            mcTH2->SetName( Form("TrueQ0VsTrueQmag{group=mecWgt,cat=%s}{group=interaction,cat=%s}", pk.mecWgt().c_str(), pk.interactionType().c_str()) );
            this->WritePlotToHistFile("TrueQ0VsTrueQmag" + pk.combination(), mcTH2, false);
            mcTH2->Draw("colz");
            for (const auto & ext : {"png", "root", "eps", "pdf"})
              c.Print( Form("%s/%s/TrueQ0VsTrueQmag_%s.%s", fOutDir.c_str(), mecWgtName.c_str(), intName.c_str(), ext) );

          }
      }

      for (const auto & plotListPair : fMCPlots)
      {
        const auto & plotName = plotListPair.first;
        const auto & plots = plotListPair.second;
//        if (plotName.find("True") == std::string::npos)
//          continue;
        if (plotName.find("Vs") != std::string::npos)
          continue;
        for (const auto & mecWgtName : MECQE_WGTS)
        {
            THStack shapeStack, sumStack;
            TLegend legend(0.7, 0.5, 0.9, 0.85);
            legend.SetFillStyle(0);
            TH1 * mcTH1 = nullptr;
            TH1 * mcTH1Area = nullptr;
            for (const auto & intClass : INT_ORDER)
            {
              PlotKey key(mecWgtName, intClass);
//              std::cout << "  trying plot: " << plotName << "/" << key << " (" << ((plots.find(key) == plots.end()) ? "NOT FOUND" : "found") << ")" << std::endl;
              mcTH1 = plots.at(key)->ToTH1(MC_NOMINAL_POT);
              mcTH1->SetLineColor(INT_COLORS.at(intClass));
              mcTH1Area = dynamic_cast<TH1*>(mcTH1->Clone());
              mcTH1Area->Scale(1./mcTH1->Integral());
              mcTH1->SetFillStyle(1001);
              mcTH1->SetFillColor(INT_COLORS.at(intClass));
              this->WritePlotToHistFile(plotName + key.combination(), mcTH1, false);
              shapeStack.Add(mcTH1Area, "hist");
              sumStack.Add(mcTH1, "hist");
              // do this the hard way to get the legend order the same as the stack instead of reversed
              legend.GetListOfPrimitives()->AddFirst(new TLegendEntry(mcTH1, intClass.c_str(), "l"));
            }
            c.Clear();
            sumStack.Draw();
            sumStack.GetXaxis()->SetTitle(mcTH1->GetXaxis()->GetTitle());
            sumStack.GetYaxis()->SetTitle("Events");
            legend.Draw();
            for (const auto & ext : {"png", "root", "eps", "pdf"})
              c.Print( Form("%s/%s/%s_MCSums.%s", fOutDir.c_str(), mecWgtName.c_str(), plotName.c_str(), ext) );

            shapeStack.Draw("nostack");
            shapeStack.GetXaxis()->SetTitle(mcTH1->GetXaxis()->GetTitle());
            shapeStack.GetYaxis()->SetTitle("Fraction of events");
            legend.Draw();
            for (const auto & ext : {"png", "root", "eps", "pdf"})
              c.Print( Form("%s/%s/%s_MCShapes.%s", fOutDir.c_str(), mecWgtName.c_str(), plotName.c_str(), ext) );
        } // for (mecWgtName)
      } // for (plotListPair)

    }

    ////////

    void WeightPlots()
    {
      const auto & wgtPlots = fMCPlots.at("weight");
      for (const auto & mecWgt : MECQE_WGTS )
      {
        TCanvas c;
        THStack stack;
        TLegend legend(0.7, 0.5, 0.9, 0.85);
        legend.SetFillStyle(0);

        std::string sameOpt = "";
        for (const auto & intType : INT_ORDER)
        {
          PlotKey key (mecWgt, intType);
    //                  std::cout << "   looking for XS weight plot: " << xsWgt << " / " << intType << std::endl;
          if (wgtPlots.find(key) == wgtPlots.end())
          {
    //                          std::cout << "     ... but couldn't find it.  skip." << std::endl;
            continue;
          }

          const auto & spectrum = wgtPlots.at(key);  // too much typing to spell out the type (ha)

          TH1 * th1 = spectrum->ToTH1(MC_NOMINAL_POT, INT_COLORS.at(intType));  // unit normalized
    //                  std::cout << th1->GetEntries() << std::endl;
          double integral = th1->Integral();
          if (integral > 0)
            th1->Scale(1./integral);

          th1->Draw(("hist" + sameOpt).c_str());
          if (sameOpt == "")
            sameOpt = " same";

          this->WritePlotToHistFile("weight" + key.combination(), th1, false);

          legend.GetListOfPrimitives()->AddFirst(new TLegendEntry(th1, key.interactionType().c_str(), "l"));
        }

        legend.Draw();

        for (const auto & ext : {"png", "root", "eps", "pdf"})
          c.Print( Form("%s/weights/%s.%s", fOutDir.c_str(), mecWgt.c_str(), ext) );
      } // for(xsWgt)
    } // WeightPlots()
};

// /////////////////////////////////////
//
//  main


void KinematicsPlots(std::string hornCurrent, std::string prodName)
{
  // too many "<image> file has been created" messages
  gErrorIgnoreLevel = kWarning;

        TH1::AddDirectory(false);
        Plotter pl(hornCurrent, prodName, "/nova/ana/users/jwolcott/2p2h_tuning/2018/plots/kinematics");
        pl.DataMCComparisons();
        pl.q0q3Panel();
        pl.WeightPlots();
        pl.TruthPlots();
}

#endif