MakeGENIEMatrix.C

Go to the documentation of this file.

#include "CAFAna/Core/Sample.h"
#include "CAFAna/Prediction/CovarianceMatrix.h"
#include "CAFAna/Prediction/PredictionInterp.h"
#include "CAFAna/Systs/XSecSystLists.h"
#include "CAFAna/Analysis/CovMxManager.h"
#include "CAFAna/Analysis/Calcs.h"
#include "NuXAna/macros/Nus18/Nus18Utilities.h"

using namespace ana;

void MakeGENIEMatrix() {

  // Run for NC FHC, ND & FD

  std::vector<covmx::Sample> samples = 
  { covmx::Sample(covmx::Analysis::kNC,
                  covmx::Polarity::kFHC,
                  covmx::Detector::kNearDet),
    covmx::Sample(covmx::Analysis::kNC,
                  covmx::Polarity::kFHC,
                  covmx::Detector::kFarDet)
  };

  // Get nominal concatenated prediction

  PredictionConcat* pred_nom = CovMxManager::GetPrediction(samples);
  pred_nom->SetBinning({ FHCNDBins(), FHCFDBins() }); // new binning schemes defined in Nus18Utilities.h
  // reeeeeally need to finalise binning schemes for each detector
  // so we can stop rebinning downstream which is super annoying

  // Get the nominal prediction as a vector (so we can draw the matrix)

  auto calc = DefaultSterileCalc(4);
  SetAngles(calc);
  TVectorD v_pred(pred_nom->NBins());
  TH1D* h_pred = pred_nom->Predict(calc).ToTH1(pred_nom->GetPOTs().back());
  for (size_t i = 0; i < pred_nom->NBins(); ++i)
    v_pred(i) = h_pred->GetBinContent(i+1);
  delete h_pred;

  // Get all the GENIE systematics (PredictionInterps for these are already made!)

  std::vector<const ISyst*> genie_systs = getAllXsecSysts_2018();

  // Instantiate an empty matrix

  CovarianceMatrix* mx = new CovarianceMatrix(pred_nom->GetPOTs().back(),
    pred_nom->GetBinning());

  // Output file!

  TFile* fout = TFile::Open("covmx.root", "recreate");
  TDirectory* hist_dir = fout->mkdir("hists");

  // Now we loop over each systematic, make a matrix and add it to the
  // originally empty matrix we defined earlier

  // Actually there are a ton of these so let's just loop over the first
  // ten to save some time

  for (size_t i = 0; i < 10; ++i) {

    const ISyst* syst = genie_systs[i];

    // Get the prediction for this systematic

    PredictionConcat* pred_syst = CovMxManager::GetPrediction(samples, syst->ShortName());
    pred_syst->SetBinning({ FHCNDBins(), FHCFDBins() });
    std::vector<std::vector<const ISyst*>> corrsysts(2, std::vector<const ISyst*>());
    pred_syst->SetSysts(new SystConcat({"nc_fhc_nd", "nc_fhc_fd"}, {syst}, corrsysts));

    // Actually make the (not empty) matrix for this syst

    CovarianceMatrix* mx_syst = new CovarianceMatrix(pred_syst, { syst }, pred_syst->GetPOTs().back());
    mx_syst->PredictCovMx(pred_nom, calc);
    TH2* h_mx = mx_syst->GetCovMxRelativeTH2(v_pred);
    std::string histname = "mxhist_" + syst->ShortName();
    hist_dir->WriteTObject(h_mx, histname.c_str());

    // Add it to our total matrix, save it, and then throw it out

    mx->AddMatrix(mx_syst);
    std::string mxname = "covmx_" + syst->ShortName();
    mx_syst->SaveTo(fout, mxname.c_str());
    delete mx_syst;

  } // for syst

  // And now save our combined matrix too!

  mx->PredictCovMx(pred_nom, calc);
  mx->SaveTo(fout, "covmx_all");

  TH2* h_mx = mx->GetCovMxRelativeWithStatsTH2(v_pred);
  hist_dir->WriteTObject(h_mx, "mxhist_all");

  delete mx;
  delete fout;

  // bye

}