plot_prediction_systs.C File Reference

#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"
#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"
#include "CAFAna/shared/Ana2018/syst_variations.h"
#include "CAFAna/shared/Ana2018/draw_plots_util.h"

Go to the source code of this file.

Functions
void	PrintLatexHeader ()
void	PrintLatexComment (TString comment="")
void	PrintLatexEnd ()
void	LtxSection (TString sec, bool newpage=true, int sublevel=0, TString comment="")
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="")

Variables
TString	mainfilename = "PlottingSysts.tex"
ofstream	mainoutput
ofstream	output
ofstream	summary
ofstream	tables
bool	isNumuAna =false
Loaders::FluxType	flux = Loaders::kFHC

Function Documentation

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

Definition at line 427 of file plot_prediction_systs.C.

References allTimeWatchdog::endl, isNumuAna, ana::MakeLatexCommandName(), and output.

{
  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 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 = ""
	)

Definition at line 59 of file plot_prediction_systs.C.

References calc, ana::CompareNDDataMCFromVector(), ana::ComparePredictionsFromVector(), om::cout, allTimeWatchdog::endl, ana::FixPlotName(), runNovaSAM::gErrorIgnoreLevel, GetCalculator(), GetExtrapolationDefs(), ana::GetSystematics2018(), isFHC, isNumuAna, ana::kCombo, ana::kEEAntiextrap, ana::kEEextrap, ana::kMEAntiextrap, ana::kMEextrap, ana::kMMAntiextrap, ana::kMMextrap, kNue1bin, kNue3bin, ana::kNuMu, ana::kNuMuNoExtrap, kNumuOpt, kNumuOther, ana::kRPARESSyst2018Test, LtxSection(), mainfilename, mainoutput, ana::MakeLatexCommandName(), ana::MakeSummaryV(), SliceLID::numuid(), output, DCS_db_parser::period, PrintLatexComment(), PrintLatexEnd(), PrintLatexHeader(), summary, systs, tables, usedumb, PandAna.Demos.tute_pid_validation::var, and copy_new_pkl_files::verbose.

{
  // 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 PrintLatexComment

(

TString

comment = ""

)

Definition at line 451 of file plot_prediction_systs.C.

References allTimeWatchdog::endl, and output.

{
  output<<"% ~~~ ~~~ ~~~ ~~~ ~~~ ~~~ ~~~ ~~~\n"<<std::endl;
  output<<"% " << comment << std::endl;
  output<<"% ~~~ ~~~ ~~~ ~~~ ~~~ ~~~ ~~~ ~~~\n"<<std::endl;
}

void PrintLatexEnd

(

)

Definition at line 458 of file plot_prediction_systs.C.

References allTimeWatchdog::endl, and mainoutput.

{
  mainoutput<<"  \\end{document}\n"<<std::endl;
}

void PrintLatexHeader

(

)

Definition at line 420 of file plot_prediction_systs.C.

References allTimeWatchdog::endl, isNumuAna, and mainoutput.

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

Variable Documentation

Loaders::FluxType flux = Loaders::kFHC

Definition at line 48 of file plot_prediction_systs.C.

Referenced by ana::nueccinc::NueCCIncCrossSectionAnalysis::AddEnhancedSample(), caf::CAFMaker::AddMCTruthToVec(), ncs::Xbeam::analyze(), ncs::Xeff::analyze(), ncs::GenieTruth::analyze(), ncs::NCAna::analyze(), genie::flux::GNuMIFluxPassThroughInfo::CalcEnuWgt(), Compare_NoExtrap(), ana::ICrossSectionAnalysis::CrossSection(), CutFlow_MC(), extract_gibuu_corrs(), FD_Data_PosComp(), FD_DataOnly_PosComp(), genie::flux::GHAKKMAtmoFlux::FillFluxHisto(), genie::flux::GFLUKAAtmoFlux::FillFluxHisto(), genie::flux::GBGLRSAtmoFlux::FillFluxHisto(), simb::MCFlux::Flux(), GenerateEventsAtFixedInitState(), genie::flux::GAtmoFlux::GenerateNext_1try(), GetCommandLineArgs(), genie::flux::GAtmoFlux::GetFlux(), ana::nueccinc::NueCCIncCrossSectionAnalysis::getFlux1D(), GetLoaders2018(), GetLoaders2020(), cvn::GetParentParticleType(), ana::BeamSystOrWeightBase::GetWeight(), ana::GeniePCASyst::GetWeight(), ana::BeamSystOrWeightBase::InitializeHistograms(), ana::GeniePCASyst::InitializeHistograms(), ana::SingleNucAnalysis::LoadFrom(), ana::NumuCC2p2hAnalysis::LoadFrom(), ana::CrossSectionSpectra::LoadFrom(), ana::CrossSectionAnalysis::LoadFrom(), ana::NumuCCIncAnalysis::LoadFrom(), ana::nueccinc::NueCCIncCrossSectionAnalysis::LoadFrom(), ana::T2KToyExperiment::LoadSpectra(), readNtuple::Loop(), main(), make_DataMCComp_numu(), make_predictions_systs(), makePrediction(), NoExtrap(), genie::flux::operator<<(), plot_xsecs(), ana::PlotAllFlux(), genie::flux::GSimpleNtpMeta::Print(), MergeCollections::MergeCollections::produce(), MergeGenCollections::MergeGenCollections::produce(), lem::MakeLibrary::produce(), evgen::GENIERockGen::produce(), evgen::GENIEGen::produce(), genie::flux::GJPARCNuFluxPassThroughInfo::Reset(), ana::NumuCCIncAnalysis::Result(), ana::NumuCC2p2hAnalysis::Result1DFluxInt(), ana::SingleNucAnalysis::Result2D(), ana::NumuCC2p2hAnalysis::Result2D(), ana::NumuCCIncAnalysis::Result2D(), runGBartolAtmoFluxDriver(), runGCylindTH1FluxDriver(), runGFlukaAtmo3DFluxDriver(), genie::flux::GAtmoFlux::SelectNeutrino(), genie::flux::GAstroFlux::SetEnergyPowLawIdx(), sk_sample_norm_abs(), caf::CAFMaker::sortNuWithIdx(), and SwapNDDataLoader().

bool isNumuAna =false

Definition at line 41 of file plot_prediction_systs.C.

TString mainfilename = "PlottingSysts.tex"

Definition at line 35 of file plot_prediction_systs.C.

ofstream mainoutput

Definition at line 36 of file plot_prediction_systs.C.

Referenced by plot_prediction_systs(), PrintLatexEnd(), and PrintLatexHeader().

ofstream output

Definition at line 37 of file plot_prediction_systs.C.

ofstream summary

Definition at line 38 of file plot_prediction_systs.C.

ofstream tables

Definition at line 39 of file plot_prediction_systs.C.