MakeCovMx.C

Go to the documentation of this file.

/// MakeCovMx.C by J. Hewes (jhewes15@fnal.gov)
// Macro to make covariance matrices

#include "CAFAna/Core/Sample.h"
#include "CAFAna/Core/ISyst.h"
#include "CAFAna/Core/HistCache.h"
#include "CAFAna/Systs/XSecSystLists.h"
#include "CAFAna/Prediction/CovarianceMatrix.h"
#include "CAFAna/Analysis/Calcs.h"
#include "NuXAna/nus/Nus19FHC/Utilities.h"

using namespace ana;

void MakeCovMx() {

  // In a utility header, we define covmx::Sample objects for the near and far
  // detector FHC NC selections. They are instantiated as kNusFHCNearDet and
  // kNusFHCFarDet – you can see they are heavily used throughout this macro
  // to keep track of the two samples.

  // Get systs and predictions
  auto nd_systs = GetSystsFromFile(kNusFHCNearDet);
  auto fd_systs = GetSystsFromFile(kNusFHCFarDet);
  IPrediction* nd_pred_systs = LoadPrediction(kNusFHCNearDet, "systs");
  IPrediction* fd_pred_systs = LoadPrediction(kNusFHCFarDet, "systs");
  //IPrediction* nd_pred_genie = LoadPrediction(kNusFHCNearDet, "genie");
  //IPrediction* fd_pred_genie = LoadPrediction(kNusFHCFarDet, "genie");

  // Create SystConcat objects to define how systematics are correlated
  // between different samples.

  // Start with a GENIE systematic which should be correlated between the
  // two detectors, for which we use a SingleSyst object. The same ISyst
  // is used across all samples, so we instantiate a SingleSyst object to
  // deal with it.
/*  SingleSyst ma_cc_res;
  ma_cc_res.syst = SystRegistry::ShortNameToSyst("MaCCRES");
  ma_cc_res.onesided = false;
  SystConcat sc_ma_cc_res({kNusFHCNearDet,kNusFHCFarDet}, {ma_cc_res}, {});
*/
  // Now we consider the Cherenkov systematic, which has a separate ISyst
  // object for each sample. Here we use the MultiSyst class to correlate
  // this systematic between the two samples. It's also a one-sided
  // systematic, so we need to tag it as such to only throw positive shifts
  // when producing the matrix.
  MultiSyst cherenkov;
  cherenkov.name = "Cherenkov";
  cherenkov.systs = {SystRegistry::ShortNameToSyst(kNusFHCNearDet.GetTag()+"_Cherenkov"),
    SystRegistry::ShortNameToSyst(kNusFHCFarDet.GetTag()+"_Cherenkov")};
  cherenkov.onesided = true;
  SystConcat sc_cherenkov({kNusFHCNearDet,kNusFHCFarDet}, {}, {cherenkov});

  // Now we do the AbsCalib systematic, which has a separate ISyst object
  // for each sample. However, we want to treat calibration as uncorrelated
  // between the two detectors, so we define a separate MultiSyst for each one
  // so they'll be thrown independently.
  MultiSyst abs_calib_nd;
  abs_calib_nd.name = "AbsCalib";
  abs_calib_nd.systs = {SystRegistry::ShortNameToSyst(kNusFHCNearDet.GetTag()+"_AbsCalib"),
    nullptr};
  abs_calib_nd.onesided = false;
  MultiSyst abs_calib_fd;
  abs_calib_fd.name = "AbsCalib";
  abs_calib_fd.systs = { nullptr,
    SystRegistry::ShortNameToSyst(kNusFHCFarDet.GetTag()+"_AbsCalib")};
  abs_calib_fd.onesided = false;
  SystConcat sc_abs_calib({ kNusFHCNearDet, kNusFHCFarDet } , {}, {abs_calib_nd,abs_calib_fd});

  // Get nominal prediction vector
  TVectorD v_nd, v_fd, v_both;
  osc::OscCalculatorSterileTrivial calc_trivial;
  TH1D* h_nd = nd_pred_systs->Predict(&calc_trivial).ToTH1(12e20);
  TH1D* h_fd = fd_pred_systs->Predict(&calc_trivial).ToTH1(12e20);
  size_t nbins_nd = h_nd->GetNbinsX(), nbins_fd = h_fd->GetNbinsX();
  v_nd.ResizeTo(nbins_nd);
  v_fd.ResizeTo(nbins_fd);
  v_both.ResizeTo(nbins_nd+nbins_fd);
  for (size_t i = 0; i < nbins_nd; ++i) {
    v_nd(i) = h_nd->GetBinContent(i+1);
    v_both(i) = h_nd->GetBinContent(i+1);
  }
  for (size_t i = 0; i < nbins_fd; ++i) {
    v_fd(i) = h_fd->GetBinContent(i+1);
    v_both(nbins_nd+i) = h_fd->GetBinContent(i+1);
  }
  HistCache::Delete(h_nd);
  HistCache::Delete(h_fd);

  auto calc = DefaultSterileCalc(4);

  // Output files
  TFile* f_out = TFile::Open("covmx.root", "recreate");
  TFile* f_plot = TFile::Open("covmx_plots.root", "recreate");

  // Produce GENIE covariance matrix
/*  CovarianceMatrix mx_ma_cc_res({nd_pred_genie,fd_pred_genie}, {12e20,12e20}, sc_ma_cc_res);
  mx_ma_cc_res.PredictCovMx({nd_pred_genie,fd_pred_genie}, {12e20,12e20}, calc);
  mx_ma_cc_res.SaveTo(f_out->mkdir("ma_cc_res"));
  TH2D* h_ma_cc_res = mx_ma_cc_res.GetCovMxRelativeTH2(v_both);
  f_plot->WriteTObject(h_ma_cc_res, "ma_cc_res");
  HistCache::Delete(h_ma_cc_res);
*/
  // Produce Cherenkov covariance matrix
  CovarianceMatrix mx_cherenkov({nd_pred_systs,fd_pred_systs}, {12e20,12e20}, sc_cherenkov);
  mx_cherenkov.PredictCovMx({nd_pred_systs,fd_pred_systs}, {12e20,12e20}, calc);
  mx_cherenkov.SaveTo(f_out->mkdir("cherenkov"));
  TH2D* h_cherenkov = mx_cherenkov.GetCovMxRelativeTH2(v_both);
  f_plot->WriteTObject(h_cherenkov, "cherenkov");
  HistCache::Delete(h_cherenkov);

  // Produce AbsCalib covariance matrix
  CovarianceMatrix mx_abs_calib({nd_pred_systs,fd_pred_systs}, {12e20,12e20}, sc_abs_calib);
  mx_abs_calib.PredictCovMx({nd_pred_systs,fd_pred_systs}, {12e20,12e20}, calc);
  mx_abs_calib.SaveTo(f_out->mkdir("abs_calib"));
  TH2D* h_abs_calib = mx_abs_calib.GetCovMxRelativeTH2(v_both);
  f_plot->WriteTObject(h_abs_calib, "abs_calib");
  HistCache::Delete(h_abs_calib);

  // Clean up
  delete nd_pred_systs;
  delete fd_pred_systs;
//  delete nd_pred_genie;
//  delete fd_pred_genie;
  delete f_out;
  delete f_plot;

}