MatrixValidation.C

Go to the documentation of this file.

/// MatrixValidation.C by J. Hewes <jhewes15@fnal.gov>

#ifdef __CINT__

#include <iostream>

void MatrixValidation(bool rhc = false) {
  throw std::runtime_error("You must run in compiled mode.");
}

#else

#include "TLine.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;

TLatex* MiscText(double x, double y, double size, TString text);
std::string GetName(int i);

void MatrixValidation(bool rhc = false) {

  // ROOT guard
  DontAddDirectory guard;

  // Preamble
  FormatCanvas();

  double nd_pot = rhc ? kAna2018SensitivityRHCNDPOT : kAna2018SensitivityFHCNDPOT;
  double fd_pot = rhc ? kAna2018RHCPOT : kAna2018FHCPOT;
  double fd_livetime = rhc ? kAna2018RHCLivetime : kAna2018FHCLivetime;

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

  osc::OscCalcSterile* calc = DefaultSterileCalc(4);
  SetAngles(calc);
  osc::OscCalcSterileTrivial* calc_trivial
    = new osc::OscCalcSterileTrivial(*calc);

  // Load and concatenate predictions
  PredictionConcat* sim = GetDefaultSimulation(rhc);

  // Load all the covariance matrices
  std::map<const ISyst*, CovarianceMatrix*> m;
  for (const ISyst* syst : GetNus18Systs(false, "base", "xsec")) {
    CovarianceMatrix* gen = LoadCovMx("covmx_fhc_xsec", "CovMx_Both_" + syst->ShortName());
    m.insert(std::pair<const ISyst*, CovarianceMatrix*>(syst, gen));
  }
  for (const ISyst* syst : GetNus18Systs(false, "base", "nonxsec")) {
    CovarianceMatrix* gen = LoadCovMx("covmx_fhc_nonxsec", "CovMx_Both_" + syst->ShortName());
    m.insert(std::pair<const ISyst*, CovarianceMatrix*>(syst, gen));
  }
  CovarianceMatrix* gen = LoadCovMx("covmx_fhc_ppfx", "CovMx_Both_PPFX");
  m.insert(std::pair<const ISyst*, CovarianceMatrix*>(nullptr, gen));
  CovarianceMatrix* matrix_combined = new CovarianceMatrix(sim, calc_trivial, {}, fd_pot, 10);
  CombineMatrices(matrix_combined, m, sim, calc);


  // Now draw these matrices
  TCanvas* c = new TCanvas();
  std::string filename = "plots/covmx/" + polarity + "_Plots.root";
  TFile* plotfile = new TFile(filename.c_str(), "recreate");

  for (auto matrix : m) {

    // Full matrix
    TH2D* h_full = matrix.second->GetFullCovMxTH2();
    int nbins = h_full->GetNbinsX();
    double last_edge = sim->GetBinEdges().back();
    double max = 13. * last_edge;
    h_full->Draw("colz");

    for (int i = 1; i < 13; ++i) {
      double pos = i * last_edge;
      TLine* lx = new TLine(0, pos, max, pos);
      lx->Draw();
      TLine* ly = new TLine(pos, 0, pos, max);
      ly->Draw();
    } 

    for (int i = 0; i < 13; ++i) {
      double pos = 0.16 + (0.72*(i/13.));
      auto l1 = MiscText(pos,0.88,0.04,GetName(i).c_str());
      l1->SetTextAngle(38);
      auto l2 = MiscText(0.12,pos,0.04,GetName(i).c_str());
      l2->SetTextAlign(31);
    }

    h_full->SetNdivisions(0,"X");
    h_full->SetNdivisions(0,"Y");

    std::string syst_name = (matrix.first == nullptr) ? "PPFX" : matrix.first->ShortName();
    std::string plot_name = "CovMx_" + polarity + "_" + syst_name + "_Full";
    std::string plot_path = "plots/covmx/" + plot_name + ".png";
    c->Write(plot_name.c_str());
    c->SaveAs(plot_path.c_str());
    c->Clear();

    // Oscillated matrix
    TH1D* nominal = sim->Predict(calc).ToTH1(fd_pot);
    nbins = nominal->GetNbinsX();
    TVectorD v_pred(nbins);
    for (int i = 0; i < nbins; ++i)
      v_pred(i) = nominal->GetBinContent(i+1);
    TH2D* h_osc = matrix.second->GetCovMxRelativeTH2(v_pred);
    h_osc->Draw("colz");
    h_osc->GetYaxis()->SetTitle("");
    h_osc->GetXaxis()->SetTitle("Deposited energy bins");
    h_osc->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");
    plot_name = "CovMx_" + polarity + "_" + syst_name + "_Oscillated";
    plot_path = "plots/covmx/" + plot_name + ".png";
    c->Write(plot_name.c_str());
    c->SaveAs(plot_path.c_str());
    c->Clear();
  }

  //
  // EFFECTIVE PULL EXTRACTION
  //

  // Get vector of all ISysts
  std::vector<const ISyst*> systs;
  std::vector<std::string> syst_types = { "xsec", "nonxsec" };
  for (std::string t : syst_types) {
    std::vector<const ISyst*> s = GetNus18Systs(rhc, "base", t);
    systs.insert(systs.end(), s.begin(), s.end());
  }

  // Loop over this vector
  for (const ISyst* syst : systs) {

    // Generate new empty matrix
    CovarianceMatrix* mx = new CovarianceMatrix(sim, calc_trivial, {}, fd_pot, 10);

    for (auto matrix : m) {

      std::string syst_name = matrix.first == nullptr ? "PPFX" : matrix.first->ShortName();

      // Drop this systematic out
      if (syst->ShortName() == syst_name) continue;

      // Add the other systematics to the new matrix
      mx->AddMatrix(matrix.second);
    }

    // Predict matrix to ensure we've reset it correctly
    mx->PredictCovMx(sim, calc);

    // Draw matrix
    auto nominal = sim->Predict(calc).ToTH1(fd_pot);
    double 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->SaveAs(Form("plots/covmx/CovMxWithout%sOscRel.png",syst->ShortName().c_str()));
    c->SaveAs(Form("plots/covmx/CovMxWithout%sOscRel.pdf",syst->ShortName().c_str()));
    c->Clear();
    delete h_osc_rel;

    // Now extract the pull manually
    SystShifts shift;
    double pull = -5;
    for (int i = 0; i < 21; ++i) {
      shift.SetShift(syst, pull);
      pull += 0.5;
    }
  }
}

//------------------------------------------------------------------------------
TLatex* MiscText(double x, double y, double size, TString text)
{
  TLatex *l = new TLatex(x, y, text);
  l->SetNDC();
  l->SetTextSize(size);
  l->SetTextColor(kBlack);
  l->Draw();
  return l;
}

//------------------------------------------------------------------------------
std::string GetName(int i) {
  if      (i==0) return std::string("#nu_{e} #rightarrow #nu_{e}");
  else if (i==1) return std::string("#bar{#nu}_{e} #rightarrow #bar{#nu}_{e}");
  else if (i==2) return std::string("#nu_{e} #rightarrow #nu_{#mu}");
  else if (i==3) return std::string("#bar{#nu}_{e} #rightarrow #bar{#nu}_{#mu}");
  else if (i==4) return std::string("#nu_{e} #rightarrow #nu_{#tau}");
  else if (i==5) return std::string("#bar{#nu}_{e} #rightarrow #bar{#nu}_{#tau}");
  else if (i==6) return std::string("#nu_{#mu} #rightarrow #nu_{e}");
  else if (i==7) return std::string("#bar{#nu}_{#mu} #rightarrow #bar{#nu}_{e}");
  else if (i==8) return std::string("#nu_{#mu} #rightarrow #nu_{#mu}");
  else if (i==9) return std::string("#bar{#nu}_{#mu} #rightarrow #bar{#nu}_{#mu}");
  else if (i==10) return std::string("#nu_{#mu} #rightarrow #nu_{#tau}");
  else if (i==11) return std::string("#bar{#nu}_{#mu} #rightarrow #bar{#nu}_{#tau}");
  else if (i==12) return std::string("NC");
  std::cerr << "What are you doing???" << std::endl;
  return "no";
}

#endif