muonid_optimization_plots.C

Go to the documentation of this file.

#pragma once

#include "CAFAna/Analysis/Exposures.h"
#include "NDAna/nuebarcc_inc/NuebarCCIncCuts.h"
#include "CAFAna/XSec/CutOptimization.h"

#include "TCanvas.h"
#include "TFile.h"
#include "TH1.h"
#include "TH1D.h"
#include "TArrow.h"
#include "TSystem.h"

using namespace ana;

void DrawArrow(TCanvas * c, double x, TString arrow_opt, int line_color = kBlack, int line_style = 1, int line_width = 3);

void muonid_optimization_plots(std::string input_file_name, std::string plot_dump)
{
  std::cout << "Loading CutOptimization found in " << input_file_name << std::endl;

  std::string results_file_name = input_file_name;
  results_file_name.insert(results_file_name.find(".root"), "_results");  
  std::cout << "Saving results to " << results_file_name << std::endl;
  std::cout << "Saving plots to " << plot_dump << std::endl;
  

  TFile * input = TFile::Open(input_file_name.c_str(), "read");
  CutOptimization * optimize = CutOptimization::LoadFrom(input, "muonid_cut_optimize")).release(;
  input->Close();

  TFile * results = new TFile(results_file_name.c_str(), "recreate");
  optimize->Optimize(CutOptimization::kIntegratedDown,
                     CutOptimization::kdSigmaOverSigma,
                     ana::kAna2019RHCPOT,
                     results->mkdir("muonid_cut_optimize"),
                     plot_dump,
                     true);

  TH1D * frac_uncert_xsec      = (TH1D*) results->Get("muonid_cut_optimize/frac_uncert_xsec"     );
  TH1D * frac_uncert_sel       = (TH1D*) results->Get("muonid_cut_optimize/frac_uncert_sel"      );
  TH1D * frac_uncert_bkgd_stat = (TH1D*) results->Get("muonid_cut_optimize/frac_uncert_bkgd_stat");
  TH1D * frac_uncert_bkgd_syst = (TH1D*) results->Get("muonid_cut_optimize/frac_uncert_bkgd_syst");
  TH1D * frac_uncert_eff       = (TH1D*) results->Get("muonid_cut_optimize/frac_uncert_eff"      );

  double maxy = 1.2;
  
  TCanvas * c = new TCanvas();
  c->SetLeftMargin(0.12);  

  frac_uncert_xsec->GetXaxis()->CenterTitle();
  frac_uncert_xsec->GetYaxis()->CenterTitle();
  frac_uncert_xsec->GetYaxis()->SetRangeUser(0, 12);
  frac_uncert_xsec->GetYaxis()->SetTitleOffset(0.8);
  frac_uncert_xsec->SetLineColor(kRed);

  frac_uncert_sel->GetXaxis()->CenterTitle();
  frac_uncert_sel->GetYaxis()->CenterTitle();
  frac_uncert_sel->GetYaxis()->SetRangeUser(0, 1.2);
  frac_uncert_sel->GetYaxis()->SetTitleOffset(0.8);
  frac_uncert_sel->SetLineColor(kRed);

  frac_uncert_bkgd_stat->GetXaxis()->CenterTitle();
  frac_uncert_bkgd_stat->GetYaxis()->CenterTitle();
  frac_uncert_bkgd_stat->GetYaxis()->SetRangeUser(0, 1.2);
  frac_uncert_bkgd_stat->GetYaxis()->SetTitleOffset(0.8);
  frac_uncert_bkgd_stat->SetLineColor(kRed);

  frac_uncert_bkgd_syst->GetXaxis()->CenterTitle();
  frac_uncert_bkgd_syst->GetYaxis()->CenterTitle();
  frac_uncert_bkgd_syst->GetYaxis()->SetRangeUser(0, 12);
  frac_uncert_bkgd_syst->GetYaxis()->SetTitleOffset(0.8);
  frac_uncert_bkgd_syst->SetLineColor(kRed);

  frac_uncert_eff->GetXaxis()->CenterTitle();
  frac_uncert_eff->GetYaxis()->CenterTitle();
  frac_uncert_eff->GetYaxis()->SetRangeUser(0, 1.2);
  frac_uncert_eff->GetYaxis()->SetTitleOffset(0.8);
  frac_uncert_eff->SetLineColor(kRed);

  frac_uncert_xsec->Draw("hist");
  c->Update();
  DrawArrow(c, 0, "<");
  c->Print((plot_dump + "/muonid_frac_uncert_xsec.pdf").c_str());

  frac_uncert_sel->Draw("hist");
  c->Update();
  DrawArrow(c, 0, "<");
  c->Print((plot_dump + "/muonid_frac_uncert_sel.pdf").c_str());

  frac_uncert_bkgd_stat->Draw("hist");
  c->Update();
  DrawArrow(c, 0, "<");
  c->Print((plot_dump + "/muonid_frac_uncert_bkgd_stat.pdf").c_str());

  frac_uncert_bkgd_syst->Draw("hist");
  c->Update();
  DrawArrow(c, 0, "<");
  c->Print((plot_dump + "/muonid_frac_uncert_bkgd_syst.pdf").c_str());

  frac_uncert_eff->Draw("hist");
  c->Update();
  DrawArrow(c, 0, "<");
  c->Print((plot_dump + "/muonid_frac_uncert_eff.pdf").c_str());

  // debug plots
  TFile * fdebug = TFile::Open((plot_dump + "/debug.root").c_str());
  auto syst_defs = optimize->GetSystDefs();
  auto mvsyst_defs = optimize->GetMVSystDefs();
  syst_defs.insert(syst_defs.end(), mvsyst_defs.begin(), mvsyst_defs.end());
  std::vector<TH1*> vabs_uncert_eff_syst;  
  std::vector<TH1*> vabs_uncert_bkgd_syst; 
  std::vector<TH1*> vfrac_uncert_eff_syst; 
  std::vector<TH1*> vfrac_uncert_bkgd_syst;
  std::vector<std::string> syst_labels;
  for(auto syst : syst_defs) {
    syst_labels.push_back(syst->GetName());
    vabs_uncert_eff_syst  .push_back((TH1*) fdebug->Get((syst->GetName() + "_abs_uncert_eff_syst"  ).c_str()));
    vabs_uncert_bkgd_syst .push_back((TH1*) fdebug->Get((syst->GetName() + "_abs_uncert_bkgd_syst" ).c_str()));
    vfrac_uncert_eff_syst .push_back((TH1*) fdebug->Get((syst->GetName() + "_frac_uncert_eff_syst" ).c_str()));
    vfrac_uncert_bkgd_syst.push_back((TH1*) fdebug->Get((syst->GetName() + "_frac_uncert_bkgd_syst").c_str()));
    vfrac_uncert_eff_syst .back()->GetYaxis()->SetRangeUser(0, 1.2);
    vfrac_uncert_bkgd_syst.back()->GetYaxis()->SetRangeUser(0, 1.2);

  }
  
  CutOptimization::PlotDebug(vabs_uncert_eff_syst  , syst_labels, "hist", "Abs. Syst. Uncert Efficiency", plot_dump + "/debug_syst_abs_uncert_eff.pdf");
  CutOptimization::PlotDebug(vabs_uncert_bkgd_syst , syst_labels, "hist", "Abs. Syst. Uncert Background", plot_dump + "/debug_syst_abs_uncert_bkgd.pdf");
  CutOptimization::PlotDebug(vfrac_uncert_eff_syst , syst_labels, "hist", "Frac. Syst. Uncert Efficiency", plot_dump + "/debug_syst_frac_uncert_eff.pdf");
  CutOptimization::PlotDebug(vfrac_uncert_bkgd_syst, syst_labels, "hist", "Frac. Syst. Uncert Background", plot_dump + "/debug_syst_frac_uncert_bkgd.pdf");

  
}

void DrawArrow(TCanvas * c, double x, TString arrow_opt, int line_color, int line_style, int line_width)
{

  bool left = arrow_opt.Contains("<");
  double offset = 0.07 * (c->GetX2() - c->GetX1());
  if(left) offset *= -1;

  double y = c->GetY1() + 0.7 * (c->GetY2() - c->GetY1());

  TLine * l1 = new TLine(x, 0, x, y);
  TArrow * arrow = new TArrow(x, y, x+offset, y, 0.04, arrow_opt.Data());
  if(left) arrow = new TArrow(x+offset, y, x, y, 0.04, arrow_opt.Data());  
  
  arrow->SetLineColor(line_color);
  l1   ->SetLineColor(line_color);

  arrow->SetLineStyle(line_style);
  l1   ->SetLineStyle(line_style);

  arrow->SetLineWidth(line_width);
  l1   ->SetLineWidth(line_width);

  arrow->Draw();
  l1->Draw();
  c->Update();
}