DrawCovMx.C

Go to the documentation of this file.

/*
 * DrawCovMx.C
 *
 *  Created on: February 15, 2018
 *      Author: J. Hewes <jhewes15@fnal.gov> 
 *
 *  Make a whole bunch o' plots :)
 */

#ifdef __CINT__

#include <iostream>
void DrawCovMx(bool rhc = false)
{
  std::cerr << "Sorry, you must run in compiled mode" << std::endl;
}

#else

#include "TH2.h"
#include "TLine.h"
#include "TLatex.h"

#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Prediction/PredictionConcat.h"
#include "CAFAna/Prediction/CovarianceMatrix.h"
#include "NuXAna/Systs/Nus18Systs.h"
#include "CAFAna/Analysis/Exposures.h"

#include "NuXAna/macros/Nus18/Nus18Utilities.h"

using namespace ana;

void DrawCovMx(bool rhc = false)
{
  DontAddDirectory guard;
  FormatCanvas();

  std::string polarity = rhc? "RHC" : "FHC";

  // Set the total POT and livetime
  double nd_pot = rhc? kAna2018SensitivityRHCNDPOT : kAna2018SensitivityFHCNDPOT;
  double fd_pot = rhc? kAna2018RHCPOT : kAna2018FHCPOT;
  double fd_livetime = rhc? kAna2018RHCLivetime : kAna2018FHCLivetime;

  // Get simulation
  PredictionConcat* sim = GetDefaultSimulation();

  // Define oscillation calculator
  osc::OscCalcSterile* calc = DefaultSterileCalc(4);
  SetAngles(calc);

  // Get matrix
  CovarianceMatrix* mx = GetFullCovMx(sim, calc, fd_pot);

  // Draw some plots
  TCanvas* c = new TCanvas();

  // First up, a bunch of covariance matrices
  auto components = GetComponents();

  TFile* plotfile = new TFile(Form("plots/covmx/%sPlots.root",polarity.c_str()),"recreate");
  
  //
  // Full covariance matrix
  //

  TH2D* h_full = mx->GetFullCovMxTH2();
  FormatFullCovMxPlot(h_full, sim);

  c->Write(Form("CovMx%sFull",polarity.c_str()));
  c->SaveAs(Form("plots/covmx/CovMx%sFull.png",polarity.c_str()));
  c->SaveAs(Form("plots/covmx/CovMx%sFull.pdf",polarity.c_str()));
  c->Clear();

  //
  // Absolute covariance matrix
  //

  std::vector<double> bin_edges = sim->GetBinEdges();

  TH2D* h_osc_abs = mx->GetCovMxAbsoluteTH2();
  h_osc_abs->Draw("colz");
  c->Write(Form("CovMx%sOscAbs",polarity.c_str()));
  c->SaveAs(Form("plots/covmx/CovMx%sOscAbs.png",polarity.c_str()));
  c->SaveAs(Form("plots/covmx/CovMx%sOscAbs.pdf",polarity.c_str()));
  c->Clear();

  //
  // Relative covariance matrix
  //

  auto nominal = sim->Predict(calc).ToTH1(fd_pot);
  int nbins = nominal->GetNbinsX();
  TVectorD v_pred(nbins);
  for (int i = 0; i < nbins; ++i)
    v_pred(i) = nominal->GetBinContent(i+1);
  TH2D* h_osc_rel = mx->GetCovMxRelativeTH2(v_pred);
  h_osc_rel->Draw("colz");
  h_osc_rel->GetYaxis()->SetTitle("");
  h_osc_rel->GetXaxis()->SetTitle("Deposited energy bins");
  h_osc_rel->GetXaxis()->CenterTitle();

  Preliminary();

  double max_e = 20;

  TLine lx(0, max_e, 2 * max_e, max_e);
  lx.Draw();
  TLine ly(max_e, 0, max_e, 2 * max_e);
  ly.Draw();
  MiscText(0.32,0.88,0.06,"ND");
  MiscText(0.68,0.88,0.06,"FD");
  MiscText(0.04,0.3,0.06,"ND");
  MiscText(0.04,0.7,0.06,"FD");
  c->Write(Form("CovMx%sRelative",polarity.c_str()));
  c->SaveAs(Form("plots/covmx/CovMx%sOscRel.png",polarity.c_str()));
  c->SaveAs(Form("plots/covmx/CovMx%sOscRel.pdf",polarity.c_str()));
  c->Clear();

  plotfile->Close();

  delete sim;
}

#endif