fiducial_optimization_plots_2d.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 "CAFAna/XSec/SystematicDef.h"

#include "TCanvas.h"
#include "TH1.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TGraph.h"
#include "TArrow.h"
#include "TSystem.h"
#include "TFile.h"
#include "TPaletteAxis.h"
#include "Utilities/rootlogon.C"
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 DrawBox(double lo, double hi, int line_color = kBlack, int line_style = 7, int line_width = 3);
void MovePalette(TH1 * hist);
void PlotSysts(std::vector<TH1*> hists, std::vector<SystematicDef*> systs, double arrow_down, double arrow_up, std::string name, double limit=0);

void fiducial_optimization_plots_2d(int axis, 
                                    std::string input_file_name = "/nova/ana/users/ddoyle/NuebarFiducialOptimization/fiducial_optimization_nocvn_all.root",
                                    std::string plot_dump       = "/nova/ana/users/ddoyle/NuebarFiducialOptimization/plots_2d_nocvn/")
{
  std::string rootdir = "fiducial_vtx_";
  double arrow_up   = 0;
  double arrow_down = 0;
  switch(axis) {
  case 0: 
    plot_dump += "x/";
    rootdir += "x";
    arrow_down = nuebarccinc::vNueCCIncFiducialMin.X();
    arrow_up   = nuebarccinc::vNueCCIncFiducialMax.X();
    break;
  case 1:
    plot_dump += "y/";
    rootdir += "y";
    arrow_down = nuebarccinc::vNueCCIncFiducialMin.Y();
    arrow_up   = nuebarccinc::vNueCCIncFiducialMax.Y();
    break;
  case 2:
    plot_dump += "z/";
    rootdir += "z";
    arrow_down = nuebarccinc::vNueCCIncFiducialMin.Z();
    arrow_up   = nuebarccinc::vNueCCIncFiducialMax.Z();
    break;    
  }
  
  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, rootdir.c_str())).release(;
  input->Close();

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

  TH2D * frac_uncert_xsec      = (TH2D*) results->Get((rootdir + "/frac_uncert_xsec"     ).c_str());
  TH2D * frac_uncert_sel       = (TH2D*) results->Get((rootdir + "/frac_uncert_sel"      ).c_str());
  TH2D * frac_uncert_bkgd_stat = (TH2D*) results->Get((rootdir + "/frac_uncert_bkgd_stat").c_str());
  TH2D * frac_uncert_bkgd_syst = (TH2D*) results->Get((rootdir + "/frac_uncert_bkgd_syst").c_str());
  TH2D * frac_uncert_eff       = (TH2D*) results->Get((rootdir + "/frac_uncert_eff"      ).c_str());

  double maxy = 1.2;
  
  TCanvas * c = new TCanvas();
  c->SetRightMargin(0.2);  

  CenterTitles(frac_uncert_xsec);
  frac_uncert_xsec->GetZaxis()->SetRangeUser(0,6);

  CenterTitles(frac_uncert_sel);
  frac_uncert_sel->GetZaxis()->SetRangeUser(0, 1.2);

  CenterTitles(frac_uncert_bkgd_stat);
  frac_uncert_bkgd_stat->GetZaxis()->SetRangeUser(0, 1.2);

  CenterTitles(frac_uncert_bkgd_syst);
  frac_uncert_bkgd_syst->GetZaxis()->SetRangeUser(0, 6);

  CenterTitles(frac_uncert_eff);
  frac_uncert_eff->GetZaxis()->SetRangeUser(0, 1.2);


  frac_uncert_xsec->Draw("colz");
  MovePalette(frac_uncert_xsec);
  c->Update();
  DrawBox(arrow_down, arrow_up);
  c->Print((plot_dump + "/" + rootdir + "_frac_uncert_xsec.pdf").c_str());

  frac_uncert_sel->Draw("colz");
  MovePalette(frac_uncert_sel);
  c->Update();
  DrawBox(arrow_down, arrow_up);
  c->Print((plot_dump + "/" + rootdir + "_frac_uncert_sel.pdf").c_str());

  frac_uncert_bkgd_stat->Draw("colz");
  MovePalette(frac_uncert_bkgd_stat);
  c->Update();
  DrawBox(arrow_down, arrow_up);
  c->Print((plot_dump + "/" + rootdir + "_frac_uncert_bkgd_stat.pdf").c_str());

  frac_uncert_bkgd_syst->Draw("colz");
  MovePalette(frac_uncert_bkgd_syst);
  c->Update();
  DrawBox(arrow_down, arrow_up);
  c->Print((plot_dump + "/" + rootdir + "_frac_uncert_bkgd_syst.pdf").c_str());

  frac_uncert_eff->Draw("colz");
  MovePalette(frac_uncert_eff);
  c->Update();
  DrawBox(arrow_down, arrow_up);
  c->Print((plot_dump + "/" + rootdir + "_frac_uncert_eff.pdf").c_str());


  // make 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;
  for(auto syst : syst_defs) {
    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()));
  }
  
  PlotSysts(vabs_uncert_eff_syst  , syst_defs, arrow_down, arrow_up, plot_dump + "/vabs_uncert_eff_syst"  );  
  PlotSysts(vabs_uncert_bkgd_syst , syst_defs, arrow_down, arrow_up, plot_dump + "/vabs_uncert_bkgd_syst" );   
  PlotSysts(vfrac_uncert_eff_syst , syst_defs, arrow_down, arrow_up, plot_dump + "/vfrac_uncert_eff_syst" , 1);  
  PlotSysts(vfrac_uncert_bkgd_syst, syst_defs, arrow_down, arrow_up, plot_dump + "/vfrac_uncert_bkgd_syst", 1);

  
}

void PlotSysts(std::vector<TH1*> hists, std::vector<SystematicDef*> systs, double arrow_down, double arrow_up, std::string name, double limit)
{
  TCanvas * c = new TCanvas(UniqueName().c_str(), "");
  c->SetRightMargin(0.2);
  for(auto ihist = 0u; ihist < hists.size(); ihist++) {
    hists[ihist]->SetTitle(systs[ihist]->GetName().c_str());
    if(limit) hists[ihist]->GetZaxis()->SetRangeUser(0, limit);
    hists[ihist]->Draw("colz");
    c->Update();
    MovePalette(hists[ihist]);
    c->Update();
    DrawBox(arrow_down, arrow_up);

    c->Print((name + "_" +systs[ihist]->GetName() + ".pdf").c_str());
  }
}

void MovePalette(TH1 * hist)
{
  auto palette = (TPaletteAxis*) hist->GetListOfFunctions()->FindObject("palette");
  palette->SetX1NDC(0.81);
  palette->SetX2NDC(0.85);
  hist->Draw("colz");
}

void DrawBox(double lo, double hi, int line_color, int line_style, int line_width)
{
  TGraph * glo = new TGraph();
  TGraph * ghi = new TGraph();
  glo->SetPoint(0, lo, -1e6);
  glo->SetPoint(1, lo,  1e6);

  ghi->SetPoint(0, -1e6, hi);
  ghi->SetPoint(1,  1e6, hi);

  glo->SetLineStyle(line_style);
  glo->SetLineWidth(line_width);
  glo->SetLineColor(line_color);
  ghi->SetLineStyle(line_style);
  ghi->SetLineWidth(line_width);
  ghi->SetLineColor(line_color);
  
  glo->Draw("same");
  ghi->Draw("same");
  
}

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();
}