plot_prediction_systs.C

Go to the documentation of this file.

// Macro to make plots comparing nominal and shifted predictions
// Run after make_nueFDprediction_shifted
// See full example in run_nueFDprediction.sh

#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/Cut.h"
#include "CAFAna/Core/HistAxis.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/Spectrum.h"
#include "3FlavorAna/Vars/NueVars.h"
#include "CAFAna/Vars/Vars.h"
#include "3FlavorAna/Prediction/PredictionSystJoint2018.h"
#include "3FlavorAna/Systs/NueSystHelper.h"  // GetSystematics2018()

#include "TFile.h"
#include "TH1.h"
#include "TCanvas.h"
#include "TGaxis.h"
#include "TLatex.h"
#include "TLegend.h"
#include "TLine.h"
#include "TSystem.h"
#include "TROOT.h"
#include "TObjectTable.h"
#include <fstream>
#include <iostream>
#include <iomanip>
#include <string>
#include <cstddef>
#include "TError.h"

// Keeping this global but redefine filename in main macro function
TString mainfilename = "PlottingSysts.tex";
ofstream mainoutput;
ofstream output;
ofstream summary;
ofstream tables;

bool isNumuAna=false; // Global variable for whether I have numu or not...

#include "CAFAna/shared/Ana2018/syst_variations.h"
#include "CAFAna/shared/Ana2018/draw_plots_util.h"

using namespace ana;

Loaders::FluxType flux = Loaders::kFHC;

void PrintLatexHeader();
void PrintLatexComment(TString comment = "");
void PrintLatexEnd();
void LtxSection (TString sec, bool newpage = true, int sublevel = 0, TString comment = "");

// For a quick test try:
//   cafe -bq plot_prediction_systs.C false FHC_FluxSyst_nue
// For debug mes try:
//   cafe -bq plot_prediction_systs.C false FHC_FluxSyst_nue allxp_nue full true
void plot_prediction_systs(
  const bool isNumu=false,
  const TString& option="FHC_All_nue",
  const TString& analysis="allxp_nue",
  const std::string& period="full",
  const bool verbose=false,
  const std::string& filename_sh="")
{
  // Set output verbosity
  // Better handling/suppressing of debug messages
  if (!verbose)
    gErrorIgnoreLevel = kWarning;

  // Messy but we need to keep isNumuAna global for now
  // TODO: make isNumuAna less global
  if (isNumu)
    isNumuAna = true;

  // Extract FHC/RHC from option
  TString HCString;
  if (option.Contains("FHC"))
  {
    isFHC = true;
    HCString = "fhc";
  }
  else
  {
    isFHC = false;
    HCString = "rhc";
  }

  const bool printtableNumuND=true;
  const bool printtableNueND=true;
  const bool printtable=true;

  // Define output name
  // Default is to use prediction file name, but may want to look at subset
  TString thisoutname (filename_sh);
  thisoutname.Remove(0,thisoutname.Last('/')+1);
  thisoutname.ReplaceAll(".root","");
  // If both option and analysis are specified construct a better name
  if (!(option == "" || analysis == ""))
  {
    thisoutname = TString(option.Data());
    thisoutname += "_" + analysis + "_" + TString(period.c_str());
  }

  // Update main filename and open file
  mainfilename = TString(thisoutname.Data());
  mainfilename += ".tex";
  mainoutput.open(mainfilename.Data());

  PrintLatexHeader();

  // Make sub-directory and files for summary and tables
  TString latexout = "latex/";
  gSystem -> MakeDirectory(latexout);
  mainoutput << "\\input{" + latexout + thisoutname + "_summary.tex}\n";
  mainoutput << "\\input{" + latexout + thisoutname + "_tables.tex}\n";
  mainoutput << "\\input{" + latexout + thisoutname + ".tex}\n";

  output.open((latexout+thisoutname+".tex").Data());
  summary.open((latexout+thisoutname+"_summary.tex").Data());
  tables.open((latexout+thisoutname+"_tables.tex").Data());

  PrintLatexComment(
    (TString)"Comparing the shifted predictions from:"+filename_sh);

  // Get info for all systs and vars
  // Define extrap type
  // Vector of extrap types to loop over
  struct xpDef
  {
    ENu2018ExtrapType xp;
    TString xpName;
    bool isExtrap;
  };
  std::vector<xpDef> xps;

  // If not doing NuMu, then fill xps_nue.
  if (!isNumuAna) {
    // No extrap
    if (analysis.Contains("nxp") || analysis.Contains("allxp"))
      xps.push_back({kNoExtrap, "pred_nxp", false});
    // Proportional decomposition
    if (analysis.Contains("xp_prop") || analysis.Contains("allxp_nue"))
      xps.push_back({kProportional, "pred_xp_prop", true});
    // Combo --> FHC only
    if (analysis.Contains("xp_combo") ||
  (analysis.Contains("allxp_nue") && isFHC))
      xps.push_back({kCombo, "pred_xp_combo", true});
  }
  // If doing NuMu, then fill xps_numu.
  else {
    // No Extrapolation
    if (analysis.Contains("nxp")    || analysis.Contains("allxp")) { xps.push_back({kNuMuNoExtrap, "pred_nxp"    , false}); }
    // Extrapolation
    if (analysis.Contains("extrap") || analysis.Contains("allxp")) { xps.push_back({kNuMu        , "pred_xp_numu", true }); }
  }

  std::cout << "Making plots for the following extrapolations: " << std::endl;
  for (const auto xp:xps)
    std::cout << "\t- " << xp.xpName << std::endl;

  // Get list of systematics
  // Pass option to determine group(s) of systematics to include in the document
  auto systs = GetSystematics2018(option);
  std::cout << "systs has size " << systs.size() << std::endl;

  // Get vector of Vars/extraps to loop over. For nue only returns Nue2018Axis now
  std::vector<ExtrapDef> xpdefs = GetExtrapolationDefs( option, period, analysis );

  for (auto xpdef:xpdefs)
    std::cout << "\t- " << xpdef.name << std::endl;
  std::cout << "" << std::endl;

  std::cout << "Making plots for the following axes (xpdefs.size is )"
    << xpdefs.size() << ": " << std::endl;
  for (auto xpdef:xpdefs)
    std::cout << "\t- " << xpdef.name << std::endl;
  std::cout << "" << std::endl;

  // Make plots and tables
  for (auto const & xp:xps)
  { // all pred types

    TString thissec = xp.isExtrap ?
      "Extrapolated prediction:":"FDMC prediction (no extrap.):";
    if(xp.xpName.Contains("prop")) thissec += " (proportional decomp.)";
    if(xp.xpName.Contains("combo")) thissec += " (combo decomp.)";
    LtxSection(thissec);

    // Define ocillation calculator
    auto calc = GetCalculator(usedumb);  // dumOsc? Default false --> global
                                         // varaible defined in
                                         // draw_plots_util

    // Make PredictionSystJoint2018 object
    bool isFakeData = false;
    if (option.Contains("fakeNDData"))
      isFakeData = true;
    std::vector<PredictionSystJoint2018*> PredictionSysts;
    if (!isNumuAna)
      PredictionSysts.push_back(
        new PredictionSystJoint2018(xp.xp, calc, systs, HCString.Data(),
                                    isFakeData, false, filename_sh, false));
    else
    {
      std::cout << "\n Looking at numu, "
        << "so first load the all quantile predictions." << std::endl;
      // First, push back the all quantile distributions.
      PredictionSysts.push_back(
        new PredictionSystJoint2018(xp.xp, calc, HCString.Data(), -1,
                                    systs, filename_sh, false));
      // Next , push back the individual quantile distributions.
      for ( int quant = 1; quant < 5; quant++ )
      {
        std::cout << "\n Now load predictions for quantile "
          << quant << std::endl;
              PredictionSysts.push_back(
          new PredictionSystJoint2018(xp.xp, calc, HCString.Data(), quant,
                                      systs, filename_sh, false));
      }
    }
    std::cout << "\n PredictionSysts has size " << PredictionSysts.size()
      << ", systs has size " << systs.size() << "\n" << std::endl;

    //auto S = systnames.begin();
    //auto shift = shifts.begin();
    //for ( ; S!=systnames.end() && shift!=shifts.end(); ++S, ++shift)
    for (const auto syst:systs)
    { // all systs

      std::cout << "Processing systematic: " << syst->ShortName() << std::endl;
      LtxSection(syst->ShortName(), true, 1,
                 "Processing systematic: " + syst->ShortName());

      // Add in something like:
      if (syst == &kRPARESSyst2018Test)
      {
        std::cout << "Skipping RPA Test..." << std::endl;
        continue;
      }

      // Loop over Vars in xpdefs
      for (auto const & var:xpdefs)
      { // all vars
        // Need to do some added switches for NuMu quantiles.
        // ApName       - Want to add "_QX" to the file names so that they are
        //                unique.
        // ThisPredSyst - For Nue simply want 0th element of PredictionSysts,
        //                but for NuMu we need to get the element which
        //                corresponds to the right quantile.
        TString ApName = "";
        PredictionSystJoint2018* ThisPredSyst = PredictionSysts[0];
        if (isNumuAna)
        {
          ApName = "_"+var.name(var.name.Index("_")+1,var.name.Sizeof());
          ApName.ReplaceAll( "Quant", "Q" );
          if      ( ApName.Contains("Q1") ) ThisPredSyst = PredictionSysts[1];
          else if ( ApName.Contains("Q2") ) ThisPredSyst = PredictionSysts[2];
          else if ( ApName.Contains("Q3") ) ThisPredSyst = PredictionSysts[3];
          else if ( ApName.Contains("Q4") ) ThisPredSyst = PredictionSysts[4];
        }

        LtxSection((syst->ShortName() + ApName), false, 2,
                   "Printing FD Plots " + MakeLatexCommandName(var.name));

        AxisType axistype = kNue1bin;
        if(var.name.Contains("Nue2018Axis")) axistype = kNue3bin;
        if(isNumuAna) axistype = kNumuOther;
        if(isNumuAna && var.name.Contains("Energy")) axistype = kNumuOpt;
        if(isNumuAna && axistype != kNumuOpt) continue;

        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        //                               // FD Plots
        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        TString AnaType = "_";
        if (!analysis.Contains("numu") && isNumuAna) AnaType = "_numu_";
        if (!analysis.Contains("nue" ) && !isNumuAna) AnaType = "_nue_";
        TString plottitle = (
          syst->ShortName() + ApName +
            (xp.isExtrap ? " FarDet Extrap (" + xp.xpName + ")":" (FarDet MC)")+
          ";" + var.name + ";" + "Events / Bin");
        TString plotname = "compare_" + xp.xpName + AnaType +
          FixPlotName(syst->ShortName()) + "_" + var.name;
        TString prefix = "nue_pred_";
        if(isNumuAna) prefix = "numu_pred_";

        bool printtable = true;

        ComparePredictionsFromVector(ThisPredSyst, syst, plottitle,
                                     plotname, printtable, axistype);

        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        //                               // ND Plots
        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        /* How do we include the ND plots using the new loading using
        PredictionSystJoint2018 concept?

        Decided with Alex that the best solution is to move the code that produces the ND plots to PredictionSystJoint2018.

        Moving towards a scenario where we have PredictionSystJoint2018 handle all the loading from file and extracting spectra.

        */
        if (!xp.isExtrap)
          continue;  // only make ND plots for extrapolated predictions

        // Define ModularExtrapComponent
        EModExtrapComps modExtrapComp;

        if (isNumuAna)
        {
          if (isFHC)
            modExtrapComp = kMMextrap;
          else
            modExtrapComp = kMMAntiextrap;
        }
        else  // Nue
        {
          if (isFHC)
            modExtrapComp = kMEextrap;
          else
            modExtrapComp = kMEAntiextrap;
        }

        TString plottitleMu = syst->ShortName() + ApName +
          "(NearDet MC); Non-Quasielastic Energy Estimator (GeV); Events / Bin";
        TString plotnameMu  = "compare_numu_ND_" +
          FixPlotName(syst->ShortName()) + "_" + var.name;
        TString numuid = "Numu";
        if(isNumuAna)
        {
          numuid = var.name;
          numuid = numuid.Remove(0, numuid.First("_")+1);
          plotnameMu = "compare_numu_ND_" + var.name + "_" +
            MakeLatexCommandName((syst->ShortName() + ApName));
        }
        LtxSection("$\\nu_{\\mu}$ ND", true, 2, "Numu Comparisons");

        CompareNDDataMCFromVector(
          ThisPredSyst, syst, modExtrapComp, plottitleMu, plotnameMu,
          numuid, kNumuOpt, printtableNumuND);
        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        if ( !printtableNueND || isNumuAna )  // continue
          continue;
        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        // Redefine ModularExtrapComponent for beam Nue background
        if (isFHC)
          modExtrapComp = kEEextrap;
        else
          modExtrapComp = kEEAntiextrap;
        TString plottitleND = syst->ShortName() + ApName +
          " (NearDet MC) "+ ";" + var.name + ";" + "Events / Bin";
        TString plotnameND  = "compare_nue_ND_"+
          FixPlotName(syst->ShortName()) + "_" + var.name;
        numuid = "Nue";
        LtxSection("$\\nu_{e}$ ND", true, 2, "Nue comparisons");
        CompareNDDataMCFromVector(
          ThisPredSyst, syst, modExtrapComp, plottitleND,
          plotnameND, numuid, kNue3bin, printtableNumuND);

      }// for var in xpdefs
    }// for S in systnames
    // Clean up PredictionSystJoint2018 objects
    for (auto predictionSyst:PredictionSysts)
    {
      delete predictionSyst;
      predictionSyst = NULL;
    }
  }// for xp in xptypes

  for (auto xp:xps)
  {
    LtxSection("Summary Tables");

    if (xp.xpName.Contains("nxp"))
    {
      LtxSection(
        "FDMC Prediction (no extrapolation) - Signal (Sorted by chisq)",
        false, 1);
      MakeSummaryV(true , "nxp", true);
      LtxSection(
        "FDMC Prediction (no extrapolation) - Background (Sorted by chisq)",
        true, 1);
      MakeSummaryV(false, "nxp", true);
    }

    if (xp.xpName.Contains("combo"))
    {
      LtxSection("Extrapolated Prediction - Signal (Sorted by chisq)", true, 1);
      MakeSummaryV(true ,"xp_combo",true);
      LtxSection("Extrapolated Prediction - Background (Sorted by chisq)",
                 true, 1);
      MakeSummaryV(false,"xp_combo",true);
    }

    if (xp.xpName.Contains("prop"))
    {
      LtxSection(
        "Extrapolated Prediction - Signal PROP (Sorted by chisq)", true, 1);
      MakeSummaryV(true ,"xp_prop",true);
      LtxSection(
        "Extrapolated Prediction - Background PROP (Sorted by chisq)", true, 1);
      MakeSummaryV(false,"xp_prop",true);
    }
    if (xp.xpName.Contains("pred_xp_numu")) {
      LtxSection  ( "Extrapolated Prediction - Signal PROP (Sorted by chisq)"    , true, 1);
      MakeSummaryV( true ,"pred_xp_numu",true);
      LtxSection  ( "Extrapolated Prediction - Background PROP (Sorted by chisq)", true, 1);
      MakeSummaryV(false ,"pred_xp_numu",true);
    }
  }

  PrintLatexEnd();

  std::cout << "Tested to HERE!" << std::endl;
  std::cout << "Finished!" << std::endl;
  return;
}

void PrintLatexHeader()
{
  if(isNumuAna) mainoutput<<"  \\input{header_numu} \n";
  else mainoutput<<"  \\input{header_nue} \n";
  mainoutput<<"  \\newpage\n"<<std::endl;
}

void LtxSection (TString sec, bool newpage, int sublevel , TString comment)
{
  if(newpage) output<<"  \\newpage\n"<<std::endl;
  output<<"\n \n% ~~~ ~~~ ~~~ ~~~ ~~~ ~~~ ~~ ~~~ ~~~ ~~~~ \n"
    << "% ~~~    " << comment
    << "% ~~~ ~~~ ~~~ ~~~ ~~~ ~~~ ~~~ ~~~ ~~~ ~~~ ~~~ \n \n"<<std::endl;

  if (isNumuAna && sublevel ==1 && !sec.Contains("Q1") )
    output<<"  \\newpage\n"<<std::endl;

  sec.ReplaceAll("_"," ");
  sec.ReplaceAll("\\nu {","\\nu_{");

  TString section = "section";
  if(sublevel ==1) section = "subsection";
  if(sublevel ==2) section = "subsubsection";
  output << "\n \\"<<section<<"{" << sec ;
  if(sublevel <=1)
    output<<"\\label{sec:" <<MakeLatexCommandName(sec) << "}";
  output << "} \n\n";

  return;
}

void PrintLatexComment(TString comment)
{
  output<<"% ~~~ ~~~ ~~~ ~~~ ~~~ ~~~ ~~~ ~~~\n"<<std::endl;
  output<<"% " << comment << std::endl;
  output<<"% ~~~ ~~~ ~~~ ~~~ ~~~ ~~~ ~~~ ~~~\n"<<std::endl;
}

void PrintLatexEnd()
{
  mainoutput<<"  \\end{document}\n"<<std::endl;
}