ppfx_syst_pca_fn.C

Go to the documentation of this file.

#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/Utilities.h"

#include "ppfx_diag_utils.h"
#include "ppfx_plot_utils.h"

#include "TCanvas.h"
#include "TH1.h"
#include "TH2.h"
#include "TFile.h"
#include "TLegend.h"
#include "TMatrixD.h"
#include "TMatrixDSymEigen.h"
#include "TMatrixDSym.h"
#include "TColor.h"
#include "TStyle.h"

#include <fstream>
#include <iostream>
#include <cmath>
#include <algorithm>


using namespace ana;

// define the flavors
std::vector<std::string> comps = {
  "nd_RHC_trueE_numu",
  "nd_RHC_trueE_nue",
  "nd_RHC_trueE_anti_numu",
  "nd_RHC_trueE_anti_nue",
  "nd_FHC_trueE_numu",
  "nd_FHC_trueE_nue",
  "nd_FHC_trueE_anti_numu",
  "nd_FHC_trueE_anti_nue",
  "fd_RHC_trueE_numu_ns",
  "fd_RHC_trueE_nue_ns",
  "fd_RHC_trueE_anti_numu_ns",
  "fd_RHC_trueE_anti_nue_ns",
  "fd_FHC_trueE_numu_ns",
  "fd_FHC_trueE_nue_ns",
  "fd_FHC_trueE_anti_numu_ns",
  "fd_FHC_trueE_anti_nue_ns",
};
std::vector<std::string> comps2 = {
  "RHC_univ_nd_trueE_numu",
  "RHC_univ_nd_trueE_nue",
  "RHC_univ_nd_trueE_anti_numu",
  "RHC_univ_nd_trueE_anti_nue",
  "FHC_univ_nd_trueE_numu",
  "FHC_univ_nd_trueE_nue",
  "FHC_univ_nd_trueE_anti_numu",
  "FHC_univ_nd_trueE_anti_nue",
  "RHC_univ_fd_trueE_numu_ns",
  "RHC_univ_fd_trueE_nue_ns",
  "RHC_univ_fd_trueE_anti_numu_ns",
  "RHC_univ_fd_trueE_anti_nue_ns",
  "FHC_univ_fd_trueE_numu_ns",
  "FHC_univ_fd_trueE_nue_ns",
  "FHC_univ_fd_trueE_anti_numu_ns",
  "FHC_univ_fd_trueE_anti_nue_ns",
};

// rebin factor for 1/32 GeV stored bin sizes
int rebin_f = 10;
int nunivs = 100;

// store flavor spectra and their universes
std::vector<TH1D*> cv_comps;
std::vector<double> cv_areas;
std::vector<std::vector<TH1D*>> univs_comps;
  
std::vector<TH1D*> cv_comps_fn;
std::vector<double> cv_areas_fn;
std::vector<std::vector<TH1D*>> univs_comps_fn;

// store shifts and the multiverse projections along each pc
std::vector<TH1D*> pc_shifts_plus;
std::vector<TH1D*> pc_shifts_minus;
// std::vector<std::vector<TH1D*>> pc_projections;

// store f/n ratios
TH1D* fn_cv;
std::vector<TH1D*> fn_shifts_ppfx;
std::vector<TH1D*> fn_shifts_plus;
std::vector<TH1D*> fn_shifts_minus;

void FillFlavorContainers(std::string fileName, std::vector<std::string> comps, std::vector<std::string> comps2) 
{
  TFile* fload = new TFile(fileName.c_str(), "READ");
  
  for(int compIdx = 0; compIdx < (int)comps.size(); ++compIdx){
    const char* comp = comps[compIdx].c_str();
    
    auto spec = Spectrum::LoadFrom(fload, TString::Format("spect_%s", comp));
    
    TH1D* temp_hist = (TH1D*)spec->ToTH1(spec->POT());
    temp_hist->Rebin(rebin_f);
    cv_areas.push_back(temp_hist->Integral());
    temp_hist->Scale(1./temp_hist->Integral());
    cv_comps.push_back(temp_hist);
  }
  
  // universes are normalized to cv areas
  for(int UnivIdx = 0; UnivIdx < nunivs; ++UnivIdx){
    std::vector<TH1D*> temp_hist_systs;
    for(int compIdx = 0; compIdx < (int)comps2.size(); ++compIdx){
      const char* comp = comps2[compIdx].c_str();
      auto spec_syst = Spectrum::LoadFrom(fload, TString::Format("spect_%s_%d", comp, UnivIdx));
      
      TH1D* hist_syst = (TH1D*)spec_syst->ToTH1(spec_syst->POT());
      hist_syst->Rebin(rebin_f);
      hist_syst->Scale(1./cv_areas[compIdx]);
      temp_hist_systs.push_back(hist_syst);
  
    }
    univs_comps.push_back(temp_hist_systs);
  }
}
 
TH1D* ConvertNormalBasis(TH1D* joint_hist_fn){
  
  std::vector<TH1D*> fn_basis = SplitHistograms(joint_hist_fn, comps.size());
  std::vector<TH1D*> near_hists(fn_basis.begin(), fn_basis.begin()+8);
  std::vector<TH1D*> fn_hists(fn_basis.begin()+8, fn_basis.end());

  std::vector<TH1D*> normal_basis;
  for(int nearIdx = 0; nearIdx < (int)near_hists.size(); nearIdx++){
    TH1D* h_normal = (TH1D*)near_hists[nearIdx]->Clone(UniqueName().c_str());
    h_normal->Scale(cv_areas_fn[nearIdx]/cv_areas[nearIdx]);
    normal_basis.push_back(h_normal);
  }
  for(int fnIdx = 0; fnIdx < (int)fn_hists.size(); fnIdx++){
    TH1D* far_h = (TH1D*)fn_hists[fnIdx]->Clone(UniqueName().c_str());
    far_h->Scale(cv_areas_fn[fnIdx+8]);
    far_h->Multiply(normal_basis[fnIdx]);
    far_h->Scale(cv_areas[fnIdx]/cv_areas[fnIdx+8]);
    normal_basis.push_back(far_h);
  }
  TH1D* ret = CombineHistograms(normal_basis);
  return ret;
}

void FillPCAContainers(TMatrixD eigenvectors, TVectorD eigenvalues, 
                       TH1D* joint_hist_cv_fn, std::vector<TH1D*> joint_hist_systs_fn,
                       TH1D* joint_hist_cv, std::vector<TH1D*> joint_hist_systs)
{
  for(int component = 0; component < eigenvalues.GetNrows(); component++){

    TVectorD x = eigenvectors[component];

    x *= sqrt(eigenvalues[component]);

    TH1D* retplus_fn = new TH1D(x);
    retplus_fn->SetName(UniqueName().c_str());
    retplus_fn->Add(joint_hist_cv_fn);
    TH1D* retplus = ConvertNormalBasis(retplus_fn);
    std::string nameplus = "pcplus"+std::to_string(component);
    retplus->SetName(nameplus.c_str());
    retplus->Divide(joint_hist_cv);
    retplus->GetYaxis()->SetRangeUser(0.7, 1.3);
    pc_shifts_plus.push_back(retplus);

    x *= -1.;

    TH1D* retminus_fn = new TH1D(x);
    retminus_fn->SetName(UniqueName().c_str());
    retminus_fn->Add(joint_hist_cv_fn);
    TH1D* retminus = ConvertNormalBasis(retminus_fn);
    std::string nameminus = "pcminus"+std::to_string(component);
    retminus->SetName(nameminus.c_str());
    retminus->Divide(joint_hist_cv);                                    // we only care about fractions
    retminus->GetYaxis()->SetRangeUser(0.7, 1.3);
    pc_shifts_minus.push_back(retminus);

  }

}

void SavePCAShifts(int ncomponents, 
                   TH1D* joint_hist_cv, std::vector<TH1D*> joint_hist_systs)
{
  TFile* saveHists = new TFile("check_cover_pc_fn_hadp_random.root", "RECREATE");

  TDirectory *cvdir = saveHists->mkdir("cv");
  cvdir->cd();
  TH1D* cv_save_hist = (TH1D*)joint_hist_cv->Clone("ppfx_cv");
  cv_save_hist->Write();
  TH1D* fn_cv_save = (TH1D*)fn_cv->Clone("fn_cv");
  fn_cv_save->Write();
  std::vector<TH1D*> hflavs = SplitHistograms(joint_hist_cv, comps.size());
  for(int comp = 0; comp < (int)comps.size(); comp++)
    hflavs[comp]->Scale(cv_areas[comp]);

  std::vector<TH1D*> near_flavs(hflavs.begin(), hflavs.begin()+8);
  std::vector<TH1D*> far_flavs(hflavs.begin()+8, hflavs.end());

  TH1D* hnear = CombineHistograms(near_flavs);
  hnear->SetName("near_cv");
  TH1D* hfar = CombineHistograms(far_flavs);
  hfar->SetName("far_cv");
  hnear->Write();
  hfar->Write();
  saveHists->cd();

  int UnivIdx=0;
  TDirectory *univdir = saveHists->mkdir("univs");
  univdir->cd();
  for(auto joint_hist_univ: joint_hist_systs) {
    std::string name_univ = "ppfx_univ"+std::to_string(UnivIdx);
    TH1D* joint_hist_syst_ratio = (TH1D*)joint_hist_univ->Clone(name_univ.c_str());
    joint_hist_syst_ratio->Write();
    UnivIdx++;
  }
  saveHists->cd();
  
  TDirectory *pcdir = saveHists->mkdir("shifts");
  pcdir->cd();
  int pcIdx = 0;
  for(auto pc_shift:pc_shifts_plus) {
    TH1D* ppfxplus = (TH1D*)pc_shift->Clone(UniqueName().c_str());
    std::string nameplus_save = "ppfx"+std::to_string(pcIdx);
    ppfxplus->SetName(nameplus_save.c_str());
    ppfxplus->Write();
    
    if(pcIdx >= ncomponents) break;
    pcIdx++;
  }
  saveHists->cd();

  TDirectory *fndir = saveHists->mkdir("fn_ratios");
  fndir->cd();
  int fnIdx = 0;
  for(auto h_fn: fn_shifts_plus) {
    TH1D* fn_ratio_save = (TH1D*) h_fn->Clone(UniqueName().c_str());
    std::string name_fn = "fn_ratios_"+std::to_string(fnIdx);
    fn_ratio_save->SetName(name_fn.c_str());
    fn_ratio_save->Write();

    if(fnIdx >= ncomponents) break;
    fnIdx++;
  }
  saveHists->cd();
  
  TDirectory *fnppfxdir = saveHists->mkdir("fn_ratios_ppfx");
  fnppfxdir->cd();
  int fnppfxIdx = 0;
  for(auto h_fn: fn_shifts_ppfx) {
    TH1D* fn_ratio_save = (TH1D*) h_fn->Clone(UniqueName().c_str());
    std::string name_fn = "fn_ratios_ppfx_"+std::to_string(fnppfxIdx);
    fn_ratio_save->SetName(name_fn.c_str());
    fn_ratio_save->Write();
    fnppfxIdx++;
  }
  saveHists->cd();
  saveHists->Close();

}

void SavePCAForCAFAna(int ncomponents, double inflate = 1.)
{
  TFile* saveHists = new TFile("ppfx_hadp_pc_shifts_fn_2018.root", "RECREATE");

  int ncomps = comps.size();
  std::vector<std::string> histnames = {
    "RHC_ND_numu",
    "RHC_ND_nue",
    "RHC_ND_anumu",
    "RHC_ND_anue",
    "FHC_ND_numu",
    "FHC_ND_nue",
    "FHC_ND_anumu",
    "FHC_ND_anue",
    "RHC_FD_numu",
    "RHC_FD_nue",
    "RHC_FD_anumu",
    "RHC_FD_anue",
    "FHC_FD_numu",
    "FHC_FD_nue",
    "FHC_FD_anumu",
    "FHC_FD_anue",
  };
  
  for(int pcIdx = 0; pcIdx < ncomponents; pcIdx++) {
    std::vector<TH1D*> hplus = SplitHistograms(pc_shifts_plus[pcIdx], ncomps); 
    std::vector<TH1D*> hminus = SplitHistograms(pc_shifts_minus[pcIdx], ncomps); 
  
    for(int hIdx = 0; hIdx < ncomps; hIdx++) {

      std::string pcIdxStr = TString::Format("%02d", pcIdx).Data();
      std::string nplus = "ppfx_hadp_pc"+pcIdxStr+"_"+histnames[hIdx]+"_plussigma";
      std::string nminus = "ppfx_hadp_pc"+pcIdxStr+"_"+histnames[hIdx]+"_minussigma";

      int nbins = hplus[hIdx]->GetNbinsX();
      TH1D* ppfxplus = new TH1D(nplus.c_str(), nplus.c_str(), nbins, 0, 10);
      TH1D* ppfxminus = new TH1D(nminus.c_str(), nminus.c_str(), nbins, 0, 10);

      for(int i = 1; i <= hplus[hIdx]->GetNbinsX(); i++){
        ppfxplus->SetBinContent(i, inflate*(hplus[hIdx]->GetBinContent(i)-1.)+1.);
        ppfxminus->SetBinContent(i, inflate*(hminus[hIdx]->GetBinContent(i)-1.)+1.);
      }
      ppfxplus->Write();
      ppfxminus->Write();
    }

  }
  saveHists->cd();
  saveHists->Close();

}

TH1D* GetFNRatio(TH1D* joint_hist)
{
  int ncomps = (int) comps.size();
  std::vector<TH1D*> hflavs = SplitHistograms(joint_hist, ncomps);
  for(int comp = 0; comp < ncomps; comp++)
    hflavs[comp]->Scale(cv_areas[comp]);

  std::vector<TH1D*> near_flavs(hflavs.begin(), hflavs.begin()+8);
  std::vector<TH1D*> far_flavs(hflavs.begin()+8, hflavs.end());
  std::vector<TH1D*> fn_flavs;
  for(int flavIdx = 0; flavIdx < (int) near_flavs.size(); flavIdx++){
    TH1D* hfn = (TH1D*)far_flavs[flavIdx]->Clone(UniqueName().c_str());
    hfn->Divide(near_flavs[flavIdx]);
    fn_flavs.push_back(hfn);
  }
  TH1D* ret = CombineHistograms(fn_flavs);

  return ret;

}

TH1D* GetTotFNRatio(TH1D* joint_hist)
{

  TH1D* hist = ConvertNormalBasis(joint_hist);
  TH1D* ret = GetFNRatio(hist);

  return ret;

}

TH1D* FNBasisHelper(TH1D* near, TH1D* far, double near_norm, double far_norm)
{
  TH1D* fn = (TH1D*)far->Clone(UniqueName().c_str());
  fn->Divide(near);
  fn->Scale(far_norm/near_norm);

  return fn;
}

void ConvertFNBasis(double nonfn_scale = 1e-6){

  for(int nearIdx = 0; nearIdx < 8; nearIdx++) {
    TH1D* near = (TH1D*)cv_comps[nearIdx]->Clone(UniqueName().c_str());
    near->Scale(nonfn_scale);
    cv_comps_fn.push_back(near);
    cv_areas_fn.push_back(cv_areas[nearIdx]/nonfn_scale);
  }
  for(int farIdx = 8; farIdx < 16; farIdx++) {
    TH1D* fn = FNBasisHelper(cv_comps[farIdx-8], cv_comps[farIdx], cv_areas[farIdx-8], cv_areas[farIdx]);
    cv_areas_fn.push_back(fn->Integral());
    fn->Scale(1./fn->Integral());
    cv_comps_fn.push_back(fn);
  }

  for(int UnivIdx = 0; UnivIdx < (int)nunivs; UnivIdx++){
    std::vector<TH1D*> temp_fn_univ;
    for(int nearIdx = 0; nearIdx < 8; nearIdx++) {
      TH1D* near_univ = (TH1D*)univs_comps[UnivIdx][nearIdx]->Clone(UniqueName().c_str());
      near_univ->Scale(nonfn_scale);
      temp_fn_univ.push_back(near_univ);
    }
    for(int farIdx = 8; farIdx < 16; farIdx++) {
      TH1D* fn_univ = FNBasisHelper(univs_comps[UnivIdx][farIdx-8], univs_comps[UnivIdx][farIdx], cv_areas[farIdx-8], cv_areas[farIdx]);
      fn_univ->Scale(1./cv_areas_fn[farIdx]);
      temp_fn_univ.push_back(fn_univ);
    }
    univs_comps_fn.push_back(temp_fn_univ);
  }
}


void ppfx_syst_pca_fn()
{
  
  // input root file      
  std::string fileName = "ppfx_smooth_flux.root";
  
  // load the spectra and normalize  
  FillFlavorContainers(fileName, comps, comps2);
  ConvertFNBasis();
  
  // flatten spectra
  TH1D* joint_hist_cv_fn = CombineHistograms(cv_comps_fn);
  std::vector <TH1D*> joint_hist_systs_fn;
  for(auto hist: univs_comps_fn) {
    TH1D* joint_hist_univ = CombineHistograms(hist);
    joint_hist_systs_fn.push_back(joint_hist_univ);
  }

  // get number of bins
  int nbins = cv_comps_fn[0]->GetNbinsX();
  int joint_nbins = joint_hist_cv_fn->GetNbinsX();

  // calculate bin-bin covariance
  const std::vector<TH1D*> ppfx(joint_hist_systs_fn.begin(), joint_hist_systs_fn.end());
  std::cout << "nunivs: " << ppfx.size();
  // std::unique_ptr<TMatrixDSym> covMx = GetCovMx(ppfx);
  std::unique_ptr<TMatrixDSym> covMx = GetCovMx(ppfx, joint_hist_cv_fn);

  // diagonalize
  TMatrixDSymEigen cov(*covMx.get());
  TMatrixD V = cov.GetEigenVectors();
  TMatrixD eigenvectors = V;
  eigenvectors.Transpose(V);
  TVectorD eigenvalues = cov.GetEigenValues();

  // Verify that the diagonalisation works as expected
  CrossCheckDiag(*covMx.get(), V, eigenvectors, eigenvalues);

  TH1D* joint_hist_cv = ConvertNormalBasis(joint_hist_cv_fn);
  std::vector<TH1D*> joint_hist_systs;
  std::for_each(joint_hist_systs_fn.begin(), joint_hist_systs_fn.end(), [&](TH1D* joint_h){
                                                                        joint_hist_systs.push_back(ConvertNormalBasis(joint_h));
               });

  // Get PCA Shifts and Projections
  FillPCAContainers(eigenvectors, eigenvalues, joint_hist_cv_fn, joint_hist_systs_fn, joint_hist_cv, joint_hist_systs);
  
  // store only these principle components
  int ncomponents = 100;

  // sum components up to get far/near ratios
  fn_cv = GetTotFNRatio(joint_hist_cv_fn);
  for(int i = 0; i < nunivs; i++)
    fn_shifts_ppfx.push_back(GetTotFNRatio(joint_hist_systs_fn[i]));

  for(int i = 0; i < ncomponents; i++){
    TH1D* fn_plus = (TH1D*)pc_shifts_plus[i]->Clone(UniqueName().c_str());
    fn_plus->Multiply(joint_hist_cv);
    fn_shifts_plus.push_back(GetFNRatio(fn_plus));

    TH1D* fn_minus = (TH1D*)pc_shifts_minus[i]->Clone(UniqueName().c_str());
    fn_minus->Multiply(joint_hist_cv);
    fn_shifts_minus.push_back(GetFNRatio(fn_minus));
  }

  SavePCAShifts(ncomponents, joint_hist_cv, joint_hist_systs);
  SavePCAForCAFAna(50, 1.25);

  // plotting!
  // store axis labels

  std::vector<std::string> axis_labels = {
    "#nu_{#mu} ND RHC", "#nu_{e} ND RHC", "#bar{#nu}_{#mu} ND RHC", "#bar{#nu}_{e} ND RHC",
    "#nu_{#mu} ND FHC", "#nu_{e} ND FHC", "#bar{#nu}_{#mu} ND FHC", "#bar{#nu}_{e} ND FHC",
    "#nu_{#mu} F/N RHC", "#nu_{e} F/N RHC", "#bar{#nu}_{#mu} F/N RHC", "#bar{#nu}_{e} F/N RHC",
    "#nu_{#mu} F/N FHC", "#nu_{e} F/N FHC", "#bar{#nu}_{#mu} F/N FHC", "#bar{#nu}_{e} F/N FHC",
  };
  std::vector<std::string> axis_labels_fn = {
    "#nu_{#mu} F/N RHC", "#nu_{e} F/N RHC", "#bar{#nu}_{#mu} F/N RHC", "#bar{#nu}_{e} F/N RHC",
    "#nu_{#mu} F/N FHC", "#nu_{e} F/N FHC", "#bar{#nu}_{#mu} F/N FHC", "#bar{#nu}_{e} F/N FHC",
  };
  std::vector<std::string> axis_labels_nom = {
    "#nu_{#mu} ND RHC", "#nu_{e} ND RHC", "#bar{#nu}_{#mu} ND RHC", "#bar{#nu}_{e} ND RHC",
    "#nu_{#mu} ND FHC", "#nu_{e} ND FHC", "#bar{#nu}_{#mu} ND FHC", "#bar{#nu}_{e} ND FHC",
    "#nu_{#mu} FD RHC", "#nu_{e} FD RHC", "#bar{#nu}_{#mu} FD RHC", "#bar{#nu}_{e} FD RHC",
    "#nu_{#mu} FD FHC", "#nu_{e} FD FHC", "#bar{#nu}_{#mu} FD FHC", "#bar{#nu}_{e} FD FHC",
  };

  // custom methods for pre-diagonalization plotting
  std::vector <plot_h> h_cv_systs;
  std::vector <std::pair<plot_h, TH1*>> h_cv_systs_ratios;

  for(auto hist: joint_hist_systs){
    plot_h univ = {hist, "hist same", kRed - 7, false, "", "l"};
    h_cv_systs.push_back(univ);
    h_cv_systs_ratios.push_back(std::make_pair(univ, joint_hist_cv));
  }

  plot_h cv = {joint_hist_cv, "hist same", kBlack, false, "", "l"};
  h_cv_systs.push_back(cv);
  h_cv_systs_ratios.push_back(std::make_pair(cv, joint_hist_cv));

  TCanvas *c1 = new TCanvas("c1", "c1", 200, 10, 800, 600);
  c1->cd();
  c1->Draw();
  plot_hists(h_cv_systs, "", false, false);
  set_labels_1D(axis_labels_nom, joint_hist_cv);
  draw_vlines(comps.size(), joint_hist_cv->GetXaxis()->GetBinLowEdge(joint_nbins+1), 0, 0.26, kBlue);
  c1->Update();
  c1->Print("ppfx_diag_plots/ppfx_hadp_spectra.pdf");

  TCanvas *c2 = new TCanvas("c2", "c2", 200, 10, 800, 600);
  c2->cd();
  c2->Draw();
  plot_ratios(h_cv_systs_ratios, "", false, 0, 2);
  set_labels_1D(axis_labels_nom, joint_hist_cv);
  draw_vlines(comps.size(), joint_hist_cv->GetXaxis()->GetBinLowEdge(joint_nbins+1), 0, 2);
  c2->Update();
  c2->Print("ppfx_diag_plots/ppfx_hadp_spectra_ratios.pdf");

  // plot the covariance matrix
  TH2D *histCov = new TH2D("histCov", "Covariance", joint_nbins, 0, joint_nbins, joint_nbins, 0, joint_nbins);
  TH2D *histCov_fn = new TH2D("histCov_fn", "Covariance", joint_nbins/2, 0, joint_nbins/2, joint_nbins/2, 0, joint_nbins/2);
  for(int xbinIdx = 1; xbinIdx <= joint_nbins; xbinIdx++){
    for(int ybinIdx = 1; ybinIdx <= joint_nbins; ybinIdx++){
      histCov->Fill(xbinIdx, ybinIdx, (*covMx)[xbinIdx-1][ybinIdx-1]);
    }
  }
  for(int xbinIdx = (joint_nbins/2)+1; xbinIdx <= joint_nbins; xbinIdx++){
    for(int ybinIdx = (joint_nbins/2)+1; ybinIdx <= joint_nbins; ybinIdx++){
      histCov_fn->Fill(xbinIdx-(joint_nbins/2), ybinIdx-(joint_nbins/2), (*covMx)[xbinIdx-1][ybinIdx-1]);
    }
  }

  gStyle->SetPalette(55);
  TCanvas *c3 = new TCanvas("c3", "c3", 200, 10, 800, 600);
  c3->cd();
  c3->Draw();
  histCov->Draw("COLZ");
  c3->SetLogz();
  c3->SetRightMargin(0.125);
  set_labels_2D(axis_labels, axis_labels, histCov);
  c3->Update();
  c3->Print("ppfx_diag_plots/ppfx_hadp_covariance.pdf");
  
  TCanvas *c3_fn = new TCanvas("c3_fn", "c3_fn", 200, 10, 800, 600);
  c3_fn->cd();
  c3_fn->Draw();
  histCov_fn->Draw("COLZ");
  c3_fn->SetLogz();
  c3_fn->SetRightMargin(0.125);
  set_labels_2D(axis_labels_fn, axis_labels_fn, histCov_fn);
  c3_fn->Update();
  c3_fn->Print("ppfx_diag_plots/ppfx_hadp_covariance_fn.pdf");

  // plot the eigenvalue profile
  TH1D* histEigen = new TH1D(eigenvalues);
  histEigen->GetXaxis()->SetTitle("Eigennumber");
  histEigen->GetYaxis()->SetTitle("Eigenvalue");
  TCanvas *c4 = new TCanvas("c4", "c4", 200, 10, 800, 600);
  c4->cd();
  c4->Draw();
  histEigen->Draw("l");
  c4->SetLogy();
  c4->Update();
  c4->Print("ppfx_diag_plots/ppfx_hadp_eigenvalues.pdf");

  int nppfx_diag_plots = 5;
  // plot the far/near ratio and shifts from the central value
  for(int component = 0; component < nppfx_diag_plots; component++){

    std::vector<plot_h> far_near_shifts = {
      {fn_cv, "hist same", kBlack, true, "CV", "l"},
      {fn_shifts_plus[component], "hist same", kBlue + 1, true, "Shift Plus", "l"},
      {fn_shifts_minus[component], "hist same", kRed, true, "Shift Minus", "l"},
    };

    std::vector<std::pair<plot_h, TH1*>> double_ratios;
    std::for_each(far_near_shifts.begin(), far_near_shifts.end(),
                    [&](plot_h h){ double_ratios.push_back(std::make_pair(h, fn_cv));
                  });

    std::vector<plot_h> pc_shifts;
    // for(int i = 0; i < (int)pc_projections[component].size(); i++)
    //   pc_shifts.push_back({pc_projections[component][i], "hist same", kGray, false, "", "l"});

    pc_shifts.push_back({pc_shifts_plus[component], "hist same", kBlue + 1, false, "", "l"});
    pc_shifts.push_back({pc_shifts_minus[component], "hist same", kRed, false, "", "l"});

    // shifts and multiverse projection ppfx_diag_plots
    TCanvas *c5 = new TCanvas("c5", "c5", 200, 10, 800, 600);
    c5->cd();
    c5->Draw();
    plot_hists(pc_shifts, "", false, false);
    plot_label(0.5, 0.93, Form("Proportion of Variance : %lg", eigenvalues[component]/eigenvalues.Sum()));
    set_labels_1D(axis_labels_nom, pc_shifts_minus[component]);
    draw_vlines(comps.size(), joint_hist_cv->GetXaxis()->GetBinLowEdge(joint_nbins+1), 0.7, 1.3);
    c5->Update();
    std::string filename1 = "ppfx_diag_plots/ppfx_hadp_shifts_pc"+std::to_string(component)+".pdf";
    c5->Print(filename1.c_str());

    TCanvas *c6 = new TCanvas("c6", "c6", 200, 10, 800, 600);
    c6->cd();
    TPad *padUp = new TPad("padUp","padUp",0.0,0.30,1.,1.00);
    padUp->Draw();
    padUp->cd();
    padUp->SetBottomMargin(0.00);
    padUp->SetBorderSize(0);

    // far/near ratio ppfx_diag_plots
    std::string fn_title = "F/N - Principle Component " + std::to_string(component);
    plot_hists(far_near_shifts, fn_title, false, false);
    // TLegend* leg = new TLegend(0.5, 0.15, 0.8, 0.35);
    // std::for_each(far_near_shifts.begin(), far_near_shifts.end(),
    //                 [&](plot_h h){ leg->AddEntry(h.hist, h.legendLabel.c_str(), h.legendOption);
    //               });
    // leg->Draw();
    padUp->Update();

    c6->cd();
    TPad *padDn = new TPad("padDn","padDn",0.0,0.0,1.,0.30);
    padDn->Draw();
    padDn->cd();
    padDn->SetBottomMargin(0.30);
    padDn->SetTopMargin(0.00);
    padDn->SetBorderSize(0);
    TH1* base = plot_ratios(double_ratios, "", false, 0.975, 1.025);
    base->GetYaxis()->SetTitle("Double Ratio");
    base->GetXaxis()->SetTitle("");
    base->GetYaxis()->SetTitleSize(base->GetYaxis()->GetTitleSize()*7./3.);
    base->GetYaxis()->SetTitleOffset(base->GetYaxis()->GetTitleOffset()*3./7.);
    base->GetYaxis()->SetLabelSize(base->GetYaxis()->GetLabelSize()*7./3.);
    base->GetXaxis()->SetTitleSize(base->GetXaxis()->GetTitleSize()*7./3.);
    base->GetXaxis()->SetTitleOffset(base->GetXaxis()->GetTitleOffset()*6./7.);
    base->GetXaxis()->SetLabelSize(base->GetXaxis()->GetLabelSize()*7./3.);
    set_labels_1D(axis_labels_fn, joint_hist_cv);
    padDn->Update();

    std::string filename2 = "ppfx_diag_plots/ppfx_hadp_fn_ratio_pc"+std::to_string(component)+".pdf";
    c6->Print(filename2.c_str());

  }
}