DataMCPair.cxx

Go to the documentation of this file.

#include "CAFAna/Analysis/DataMCPair.h"


#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/Var.h"
#include "CAFAna/Core/SpectrumLoader.h"

#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Analysis/Style.h"

#include "TDirectory.h"
#include "TH1.h"
#include "TH2.h"
#include "TKey.h"
#include "TLatex.h"
#include "TLegend.h"
#include "TList.h"
#include "TObjString.h"
#include "TParameter.h"

#include <iostream>

namespace ana
{

  // Replace all instances of from with to in str
  std::string replaceAll(const std::string& str,
                         const std::string& from,
                         const std::string& to) {
      std::string ret(str);
      if(from.empty())
          return ret;
      size_t start_pos = 0;
      while((start_pos = ret.find(from, start_pos)) != std::string::npos) {
          ret.replace(start_pos, from.length(), to);
          start_pos += to.length(); // In case 'to' contains 'from', like replacing 'x' with 'yx'
      }
      return ret;
  }


  TString TableNum(float num)
  {
    return TString::Format("%.02lf", num);
  }

  //-------------------------------------------------------------------------
  DataMCPair::DataMCPair(Selection sel, TangibleAxis tanAxis,
                         SpectrumLoaderBase& loaderData,
                         SpectrumLoaderBase& loaderMC,
                         const Var & wei,
                         const std::vector<const ISyst*>& systs,
                         const std::map<std::string, SpectrumComponents::Component>& mcCategories)
  :
    fSel(sel),
    fData(loaderData, tanAxis.fObj, fSel.fObj, kNoShift),
    fMC(fSel.fObj, tanAxis.fObj, loaderMC, wei, mcCategories),
    fShortName(tanAxis.fShortName + "_" + fSel.fShortName),
    fSysts(systs)
  {
    fUps.reserve(systs.size());
    fDowns.reserve(systs.size());
    for(const auto& syst:systs){
      fUps.emplace_back(loaderMC, tanAxis.fObj, fSel.fObj, SystShifts(syst,+1), wei);
      fDowns.emplace_back(loaderMC, tanAxis.fObj, fSel.fObj, SystShifts(syst,-1), wei);
    }
  }

  //-------------------------------------------------------------------------

  DataMCPair::DataMCPair(Cut cut, HistAxis axis,
                         SpectrumLoaderBase& loaderData,
                         SpectrumLoaderBase& loaderMC,
                         const Var & wei,
                         const std::vector<const ISyst*>& systs,
                         const std::map<std::string, SpectrumComponents::Component>& mcCategories)
   : DataMCPair(Selection(cut, "", ""), TangibleAxis(axis, "", ""),
               loaderData, loaderMC, wei, systs, mcCategories)
  {
  }


  //-------------------------------------------------------------------------
  void DataMCPair::DrawMCSyst(const int iSyst, EBinType bintype) const
  {
 
    if(iSyst < 0)
      PlotWithSystErrorBand(fMC.Tot(), fUps, fDowns, fData.POT(),
                kFullPredColor, kPredErrColor, 1.5, true, bintype);
    else
    {
      assert(iSyst < (int)fUps.size());
      PlotWithSystErrorBand(fMC.Tot(), {fUps[iSyst]}, {fDowns[iSyst]}, fData.POT(),
                kFullPredColor, kPredErrColor, 1.5, true, bintype);
    }
  }

  //-------------------------------------------------------------------------
  void DataMCPair::OverlayDataMCSyst(const int iSyst, bool drawMCBreakdown, EBinType bintype, bool drawLegend) const
  {
    DrawMCSyst(iSyst, bintype);
    if (drawMCBreakdown)
      DrawMCComponents(true, bintype);
    DrawData(kBlack, bintype);
    if (drawLegend)
      DrawLegend();
  }

  //-------------------------------------------------------------------------
  void DataMCPair::DrawMCNormSyst(const int iSyst, EBinType bintype) const
  {
    // New vectors for area normalized systs
    std::vector<Spectrum> normUps;
    std::vector<Spectrum> normDowns;

    if (iSyst >= 0)
    {
      // We'll have one entry if there's just one syst to plot.
      normUps.reserve(1);
      normDowns.reserve(1);
    }
    else
    {
      // Otherwise we need all systs.
      normUps.reserve(fUps.size());
      normDowns.reserve(fUps.size());
    }
    for(size_t i = 0; i < fUps.size(); ++i)
    {
      // bail if we want a particular syst and it's not the one we want.
      if(iSyst >= 0 && (int)i != iSyst) continue;
      normUps.emplace_back(fUps[i]);
      Spectrum& norm = normUps.back();
      norm.OverridePOT(norm.POT() * norm.Integral(1) / fMC.Tot().Integral(1));
    }
    for(size_t i = 0; i < fDowns.size(); ++i)
    {
      // bail if we want a particular syst and it's not the one we want.
      if(iSyst >= 0 && (int)i != iSyst) continue;
      normDowns.emplace_back(fDowns[i]);
      Spectrum& norm = normDowns.back();
      norm.OverridePOT(norm.POT() * norm.Integral(1) / fMC.Tot().Integral(1));
    }
    PlotWithSystErrorBand(fMC.Tot(), normUps, normDowns, fData.POT(),
              kFullPredColor, kPredErrColor, 1.5,true, bintype);
  }

  //-------------------------------------------------------------------------
  void DataMCPair::OverlayDataMCSystNorm(const int iSyst, bool drawMCBreakdown, EBinType bintype, bool drawLegend) const
  {
    DrawMCNormSyst(iSyst, bintype);
    if (drawMCBreakdown)
      DrawMCComponents(true, bintype);
    DrawData(kBlack, bintype);
    if (drawLegend)
      DrawLegend();
  }

  //-------------------------------------------------------------------------
  TLegend* DataMCPair::DrawLegend(const int dataColor, bool includeComponents, double drawThreshold) const
  {
    TLegend* leg;
    if (includeComponents)
            // needs to be a bit smaller if components included
            leg = AutoPlaceLegend(0.25, 0.18 + 0.04*fMC.ComponentDescrs().size());
    else
            leg = AutoPlaceLegend(0.36, 0.18, 0.67);
    //TLegend* leg = new TLegend(0.6, 0.7, 0.85, 0.85);
    TH1F* colMC = new TH1F();
    colMC->SetLineColor(kFullPredColor);
    TH1F* colMCBkg = new TH1F();
    colMCBkg->SetLineColor(kTotalBkgColor);
    TH1F* colData = new TH1F();
    colData->SetLineColor(dataColor);
    colData->SetMarkerColor(dataColor);

    leg->AddEntry((TObject*)colData, "Data", "ple");
    leg->AddEntry((TObject*)colMC, (includeComponents ? "Total Simulation" : "Simulation"), "l");

    if (includeComponents)
      fMC.DrawLegend(drawThreshold, leg);


    leg->SetFillColor(0);
    leg->SetFillStyle(0);

    for(const auto& obj:*leg->GetListOfPrimitives())
    {
      if(obj->InheritsFrom("TAttFill")){
        ((TAttFill*)obj)->SetFillStyle(0);
        //((TAttFill*)obj)->SetFillColor(0);
      }

    }
    leg->Draw();
    return leg;
  }

  //-------------------------------------------------------------------------
  void DataMCPair::AddExposure() const
  {
    std::stringstream expo;
    float pot = fData.POT();
    expo.precision(2);
    expo  << "ND, "
          << TString::Format("%.02lf#times10^{20}", pot/1e20) << " POT";
    TLatex text;
    text.SetTextAlign(12);

    text.DrawLatexNDC(0.2, 0.95, expo.str().c_str());
  }

  //-------------------------------------------------------------------------
  void DataMCPair::DrawData(const int color, EBinType bintype) const
  {
    TH1D* hData = fData.ToTH1(fData.POT(), color, kFullCircle, kPOT, bintype);

    hData->Draw("SAME");
  }

  //-------------------------------------------------------------------------
  void DataMCPair::DrawMCComponents(bool stacked, EBinType bintype) const
  {
    fMC.DrawComponents(stacked, bintype, fData.POT());
  }

  //-------------------------------------------------------------------------
  float DataMCPair::Purity(const std::set<std::string>& signalCatNames) const
  {
    return fMC.Purity(signalCatNames);
  }

  //-------------------------------------------------------------------------
  void DataMCPair::CreateSystTable(const std::string& fname)const
  {

    std::ofstream table(fname);
    table <<
    "\\documentclass[12pt]{article}\n"
    "\\usepackage{longtable}\n"
    "\\begin{document}\n"
    << fShortName <<
    "\\begin{longtable}\n"
    "{|p{0.4\\textwidth}|p{0.3\\textwidth}|p{0.3\\textwidth}|}\n"
    "\\hline\n"
    "Parameter & $+1 \\sigma$ Mean & $-1 \\sigma$ Mean  \\\\\n"
    "\\hline \\hline\n"
    "\\endfirsthead\n"
    "\\multicolumn{3}{l}\n"
    "{{\\bfseries\\tablename\\ \\thetable{} -- "
    "continued from previous page}}\\\\\n"
    "\\hline\n"
    "Parameter & $+1 \\sigma$ Mean & $-1 \\sigma$ Mean  \\\\\n"
    "\\hline \\hline\n"
    "\\endhead\n"
    "\\hline \\hline \\multicolumn{3}{|r|}{{Continued on next page}}"
    " \\\\ \\hline\n"
    "\\endfoot\n"
    "\\hline\n"
    "\\endlastfoot\n"
    "\\hline\n"
    << std::endl;

    float meanNom = fMC.Tot().Mean();
    for(size_t iSyst = 0; iSyst < fSysts.size(); ++ iSyst)
    {

      float meanUp = fUps[iSyst].Mean();
      float meanDown = fDowns[iSyst].Mean();
      float percUp = 100 * (meanUp - meanNom)/meanNom;
      float percDown = 100 * (meanDown - meanNom)/meanNom;

      table << replaceAll(fSysts[iSyst]->LatexName(), "_", "\\_") << " & "
            << TableNum(meanUp) << " (" << TableNum(percUp) << "\\%) &"
            << TableNum(meanDown) << " (" << TableNum(percDown) << "\\%) "
            << " \\\\ \\hline" << std::endl;
    }
    table << "\\end{longtable}\n\\end{document}" << std::endl;
  }

  //-------------------------------------------------------------------------
  void DataMCPair::SaveTo(TDirectory* dir, const std::string& name) const
  {
    TDirectory* tmp = gDirectory;

    dir = dir->mkdir(name.c_str()); // switch to subdir
    dir->cd();

    fData.SaveTo(dir, "data");

    fMC.SaveTo(dir, "MC");

    TDirectory * upsDir = dir->mkdir("upShifts");
    for (std::size_t idx = 0; idx < fUps.size(); idx++)
    {
      fUps[idx].SaveTo(upsDir, std::to_string(idx));
    }
    TDirectory * downsDir = dir->mkdir("downShifts");
    for (std::size_t idx = 0; idx < fDowns.size(); idx++)
    {
      fDowns[idx].SaveTo(downsDir, std::to_string(idx));
    }

    dir->cd();
    TObjString shortName(fShortName.c_str());
    shortName.Write("shortName");

    dir->Write();
    delete dir;

    tmp->cd();
  }

  //-------------------------------------------------------------------------
  std::unique_ptr<DataMCPair> DataMCPair::LoadFrom(TDirectory* dir, const std::string& name)
  {
    dir = dir->GetDirectory(name.c_str()); // switch to subdir
    assert(dir);

    Spectrum data = *Spectrum::LoadFrom(dir, "data");
    SpectrumComponents mc = *SpectrumComponents::LoadFrom(dir, "MC");

    std::map<std::string, std::vector<Spectrum>> specs;
    for (const auto & shift : {"up", "down"})
    {
      TDirectory * shiftsDir = dynamic_cast<TDirectory*>(dir->Get(Form("%sShifts", shift)));
      // maybe there weren't any shifts, and a hadd didn't create the empty dir?
      if (!shiftsDir)
        continue;

      for ( TObject * obj : *(shiftsDir->GetListOfKeys()) )
      {
        TKey * key = dynamic_cast<TKey*>(obj);
  
        // pray that they're read out in the same order they went in
        specs[shift].emplace_back( *(Spectrum::LoadFrom(shiftsDir, key->GetName()).release()) );
      }
    }
    auto ret = std::make_unique<DataMCPair>(
        std::move(data),
        std::move(mc),
        std::move(specs["up"]),
        std::move(specs["down"]),
        dynamic_cast<TObjString*>(dir->Get("shortName"))->String().Data()
    );

    delete dir;

    return std::move(ret);
  }

  //-------------------------------------------------------------------------

  void DataMCPair::SetComponentColor(const std::string &name, int color)
  {
    assert (fMC.ComponentDescrs().find(name) != fMC.ComponentDescrs().end());

    fMC.ComponentDescrs().at(name).color = color;
  }

  //-------------------------------------------------------------------------

  void DataMCPair::SetComponentBlurb(const std::string &name, const std::string &blurb)
  {
    assert (fMC.ComponentDescrs().find(name) != fMC.ComponentDescrs().end());

    fMC.ComponentDescrs().at(name).blurb = blurb;

  }


} // namespace