DrawSystShifts.C

Go to the documentation of this file.

// DrawSystShifts by J. Hewes (jhewes15@fnal.gov)
/// Draw systematic shifts in each detector bin and compare to covariance matrix

#include "CAFAna/Core/ISyst.h"
#include "CAFAna/Systs/XSecSystLists.h"
#include "CAFAna/Prediction/CovarianceMatrix.h"
#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Exposures.h"

#include "NuXAna/macros/Nus20/Utilities.h"

#include "TSystem.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TRandom3.h"
#include "TText.h"

#include <fstream>
#include <algorithm>

using namespace ana;

float Gaus(float mean, float var, float x) {
  float norm = 1.0/sqrt(2.0*M_PI*var);
  return norm * exp(-1*pow(x-mean,2)/(2*var));
}

void DrawBins(IPrediction* pred, const ISyst* syst, TMatrixD mx,
  TVectorD vec, size_t nbins, std::ofstream& f, std::string det);

void DrawPulls(IPrediction* pred, const ISyst* syst, size_t nbins,
  std::ofstream& f, std::string det);

size_t nPoints = 1001;
std::vector<std::string> one_d = { "AbsCalibShape", "Cherenkov", "XSecTune" };

void DrawSystShifts() {

  DontAddDirectory guard;

  // Define exposures
  double potND = kAna2018SensitivityFHCNDPOT;
  double potFD = kAna2018FHCPOT;

  // Set up latex output
  gSystem->MakeDirectory("tex_pulls");
  gSystem->MakeDirectory("tex_pulls/plots");
  gSystem->CopyFile(Form("%s/nus/Nus20/Plots/header_nus.tex",
    FindCAFAnaDir().c_str()), "tex_pulls/header_nus.tex");

  // Write header
  std::ofstream maintex;
  maintex.open("tex_pulls/main.tex");
  maintex << "\n\\input{header_nus.tex}\n\n";
  maintex << "\\newpage\n\n";

  // Load predictions
  IPrediction* predNDSysts = LoadPrediction(kNusFHCNearDet, "systs", kUseCVMFS);
  IPrediction* predFDSysts = LoadPrediction(kNusFHCFarDet, "systs", kUseCVMFS);
  IPrediction* predNDGENIE = LoadPrediction(kNusFHCNearDet, "genie", kUseCVMFS);
  IPrediction* predFDGENIE = LoadPrediction(kNusFHCFarDet, "genie", kUseCVMFS);
  auto calc = DefaultSterileCalc(4);
  SetNus20Params(calc);

  // Spectrum vectors
  TVectorD vND(kNCNDBinning.NBins());
  TVectorD vFD(kNCFDBinning.NBins());
  TH1D* hND = predNDSysts->Predict(calc).ToTH1(potND);
  TH1D* hFD = predFDSysts->Predict(calc).ToTH1(potND);
  for (size_t i = 0; i < size_t(kNCNDBinning.NBins()); ++i)
    vND(i) = hND->GetBinContent(i+1);
  for (size_t i = 0; i < size_t(kNCFDBinning.NBins()); ++i)
    vFD(i) = hFD->GetBinContent(i+1);
  HistCache::Delete(hND);
  HistCache::Delete(hFD);

  // Load covariance matrices
  std::map<std::string, CovarianceMatrix*> mxND;
  for (auto itMx : GetMatrices({kNusFHCNearDet}, kUseCVMFS))
    mxND[itMx->GetName()] = itMx;
  std::map<std::string, CovarianceMatrix*> mxFD;
  for (auto itMx : GetMatrices({kNusFHCFarDet}, kUseCVMFS))
    mxFD[itMx->GetName()] = itMx;

  // Near detector systematics
  maintex << "\\section{Near Detector}\n\n";
  for (auto syst : LoadSysts(kNusFHCNearDet, kUseCVMFS)) {
    maintex << "\\subsection{" << GetSystBaseName(syst->ShortName()) << "}\n\n";
    TMatrixD mx(mxND[GetSystBaseName(syst->ShortName())]->GetCovMxRelative(vND));
    DrawBins(predNDSysts, potND, syst, mx, vND, kNCNDBinning.NBins(), maintex, "nd");
    DrawPulls(predNDSysts, potND, syst, kNCNDBinning.NBins(), maintex, "nd");
  }
  for (auto syst : getAllXsecSysts_2018()) {
    std::string systname = syst->ShortName();
    std::replace(systname.begin(), systname.end(), '_', ' ');
    maintex << "\\subsection{" << systname << "}\n\n";
    TMatrixD mx(mxND[syst->ShortName()]->GetCovMxRelative(vND));
    DrawBins(predNDGENIE, potND, syst, mx, vND, kNCNDBinning.NBins(), maintex, "nd");
    DrawPulls(predNDGENIE, potND, syst, kNCNDBinning.NBins(), maintex, "nd");
  }

  // Far detector systematics
  maintex << "\\section{Far Detector}\n\n";
  for (auto syst : LoadSysts(kNusFHCFarDet, kUseCVMFS)) {
    maintex << "\\subsection{" << GetSystBaseName(syst->ShortName()) << "}\n\n";
    TMatrixD mx(mxFD[GetSystBaseName(syst->ShortName())]->GetCovMxRelative(vFD));
    DrawBins(predFDSysts, potFD, syst, mx, vFD, kNCFDBinning.NBins(), maintex, "fd");
    DrawPulls(predFDSysts, potFD, syst, kNCFDBinning.NBins(), maintex, "fd");
  }
  for (auto syst : getAllXsecSysts_2018()) {
    std::string systname = syst->ShortName();
    std::replace(systname.begin(), systname.end(), '_', ' ');
    maintex << "\\subsection{" << systname << "}\n\n";
    TMatrixD mx(mxFD[syst->ShortName()]->GetCovMxRelative(vFD));
    DrawBins(predFDGENIE, potFD, syst, mx, vFD, kNCFDBinning.NBins(), maintex, "fd");
    DrawPulls(predFDGENIE, potFD, syst, kNCFDBinning.NBins(), maintex, "fd");
  }

  maintex << "\\end{document}\n";
  maintex.close();

}

void DrawBins(IPrediction* pred, double pot, const ISyst* syst, TMatrixD mx,
  TVectorD vec, size_t nbins, std::ofstream& f, std::string det) {

  std::cout << "Processing syst " << syst->ShortName() << std::endl;
  f << "\\begin{figure}[h!]\n  \\centering\n";
  auto calc = DefaultSterileCalc(4);
  SetNus20Params(calc);

  // Make output dirs
  std::string systname = syst->ShortName().substr(0,7) == "nus_fhc" ?
    GetSystBaseName(syst->ShortName()) : syst->ShortName();

  std::vector<TH1D*> bins(nbins);
  bool allNoEffect = true;
  bool allOneDirectional = true;
  std::vector<bool> oneDirectional(nbins, false);
  for (size_t i = 0; i < nbins; ++i)
    bins[i] = new TH1D("", "", 60, 1-(3*sqrt(mx(i,i))), 1+(3*sqrt(mx(i,i))));

  TH1D* hNom = pred->Predict(calc).ToTH1();

  // Loop over bins
  TRandom3 rnd(0);
  for (size_t i = 0; i <= 1e5; ++i) {
    SystShifts shifts;
    if (std::any_of(one_d.begin(), one_d.end(), [systname](std::string name)
      { return name == systname; })) shifts.SetShift(syst, fabs(rnd.Gaus()));
    else shifts.SetShift(syst, rnd.Gaus());
    TH1D* h = pred->PredictSyst(calc, shifts).ToTH1(pot);
    for (size_t j = 0; j < nbins; ++j)
      bins[j]->Fill(h->GetBinContent(j+1)/vec(j));
    HistCache::Delete(h);
  }

  // Here we want to classify each systematic in each bin.
  // Here's how this is going to go. First we want to classify
  // as one-directional, or not. This should be super easy.
  // In each bin, we add up bins below and above centre, and if
  // either integral is zero, it's one-directional. That's
  // all there is to it.
  for (size_t i = 0; i < nbins; ++i) {
    if (bins[i]->Integral(1, 29) != 0 or bins[i]->Integral(32, 60) != 0)
      allNoEffect = false;
    if (bins[i]->Integral(1, 29) == 0 or bins[i]->Integral(32, 60) == 0)
      oneDirectional[i] = true;
    else allOneDirectional = false;
  }

  // Calculate area fraction. For this we literally just loop over a bunch
  // of slices in x, and take the minimum height at each one. then we add
  // up a little area slice.
  std::vector<float> nonGaussianScore(nbins, 0);
  for (size_t i = 0; i < nbins; ++i) {
    float xrange = 1 - bins[i]->GetXaxis()->GetBinLowEdge(1);
    float norm = 0.1*bins[i]->Integral()*sqrt(mx(i,i));
    for (int j = -999; j <= +999; ++j) {
      float x = (0.001 * float(j) * xrange)+1;
      float yHist = bins[i]->GetBinContent(bins[i]->GetXaxis()->FindBin(x));
      float yGaus = norm * Gaus(1, mx(i,i), x);
      if (oneDirectional[i]) yGaus *= 2;
      float yMin = yHist < yGaus ? yHist : yGaus;
      nonGaussianScore[i] += (yMin * 0.1 * xrange) / norm;
    }
  }

  // Save histograms to file
  TCanvas c("", "", 900, 1600);
  std::vector<TPad*> pads;
  std::vector<TGraph*> graphs;
  pads.push_back(new TPad("classification", "", 0, 0.95, 1, 0.999));
  pads.back()->cd();
  std::string classification;
  if (allNoEffect) classification = "No effect";
  else if (allOneDirectional) classification = "One-directional";
  else classification = "Two-directional";
  float meanScore = 0.;
  for (size_t i = 0; i < nbins; ++i) meanScore += nonGaussianScore[i] / float(nbins);
  TText* t = new TText(0.5, 0.5, Form("%s, score = %.1f%%", classification.c_str(), meanScore));
  t->SetTextSize(0.6);
  t->SetTextAlign(22);
  t->Draw();
  for (size_t i = 0; i < nbins; ++i) {
    float x1 = 0.25*(float(i%4));
    float y1 = 0.95-(0.19*(1.+std::floor(0.25*float(i))));
    pads.push_back(new TPad(Form("bin%zu", i), "", x1, y1, x1+0.25, y1+0.18));
    pads.back()->cd();
    std::vector<float> x(nPoints);
    std::vector<float> y(nPoints);
    float binWidth = 6.0*sqrt(mx(i,i))/float(nPoints-1);
    float norm = 0.1*bins[i]->Integral()*sqrt(mx(i,i));
    for (size_t j = 0; j < nPoints; ++j) {
      x[j] = 1 + (-3.0*sqrt(mx(i,i))) + (float(j)*binWidth);
      y[j] = norm * Gaus(1,mx(i,i),x[j]);
      if (oneDirectional[i]) y[j] *= 2;
    }
    graphs.push_back(new TGraph(nPoints, &x[0], &y[0]));
    bins[i]->Draw("hist");
    bins[i]->SetTitle(Form("Bin %zu - %.2f < E (GeV) < %.2f", i+1,
      hNom->GetXaxis()->GetBinLowEdge(i+1), hNom->GetXaxis()->GetBinUpEdge(i+1)));
    graphs.back()->Draw("c same");
    graphs.back()->SetLineColor(kRed-7);
    graphs.back()->SetLineWidth(3);
    TText* tmp = new TText(0.85, 0.85, Form("%.1f%%", nonGaussianScore[i]));
    tmp->SetTextAlign(33);
    tmp->SetTextSize(0.1);
    tmp->SetNDC(true);
    tmp->Draw();
    // pads.back()->Update();
  }

  c.cd();
  for (auto pad : pads) pad->Draw();

  c.SaveAs(Form("tex_pulls/plots/gaus_%s_%s.pdf",
    det.c_str(), systname.c_str()));
  f << "  \\includegraphics[height=0.8\\textheight]{"
    << Form("plots/gaus_%s_%s.pdf", det.c_str(), systname.c_str())
    << "}\n\\end{figure}\n\n\\clearpage\n\n";

  HistCache::Delete(hNom);
  delete calc;
}

void DrawPulls(IPrediction* pred, double pot, const ISyst* syst, size_t nbins,
  std::ofstream& f, std::string det) {

  f << "\\begin{figure}[h!]\n  \\centering\n";
  auto calc = DefaultSterileCalc(4);
  SetNus20Params(calc);

  // Make output dirs
  std::string systname = syst->ShortName().substr(0,7) == "nus_fhc" ?
    GetSystBaseName(syst->ShortName()) : syst->ShortName();

  std::vector<float> x(nPoints);
  std::vector<std::vector<float>> y(nbins, std::vector<float>(nPoints));

  TH1D* hNom = pred->Predict(calc).ToTH1(pot);

  for (size_t i = 0; i < nPoints; ++i) {
    SystShifts shifts;
    x[i] = (6./float(nPoints-1))*(float(i)-(0.5*(nPoints-1)));
    shifts.SetShift(syst, x[i]);
    TH1D* h = pred->PredictSyst(calc, shifts).ToTH1(pot);
    for (size_t j = 0; j < nbins; ++j)
      y[j][i] = h->GetBinContent(j+1);
    HistCache::Delete(h);
  }

  // Save histograms to file
  TCanvas c("", "", 900, 1600);
  std::vector<TPad*> pads;
  std::vector<TGraph*> graphs;
  for (size_t i = 0; i < nbins; ++i) {
    float x1 = 0.25*(float(i%4));
    float y1 = 1-(0.2*(1.+std::floor(0.25*float(i))));
    pads.push_back(new TPad(Form("bin%zu", i), "", x1, y1, x1+0.25, y1+0.1999));
    pads.back()->cd();   
    graphs.push_back(new TGraph(nPoints, &x[0], &y[i][0]));
    graphs.back()->Draw("al");
    graphs.back()->SetTitle(Form("Bin %zu - %.2f < E (GeV) < %.2f", i+1,
      hNom->GetXaxis()->GetBinLowEdge(i+1), hNom->GetXaxis()->GetBinUpEdge(i+1)));
    graphs.back()->SetLineWidth(3);
  }

  c.cd();
  for (auto pad : pads) pad->Draw();

  c.SaveAs(Form("tex_pulls/plots/pull_%s_%s.pdf",
    det.c_str(), systname.c_str()));
  f << "  \\includegraphics[height=0.8\\textheight]{"
    << Form("plots/pull_%s_%s.pdf", det.c_str(), systname.c_str())
    << "}\n\\end{figure}\n\n\\clearpage\n\n";

  delete calc;
}