genie_syst_pca.C

Go to the documentation of this file.

#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Analysis/Exposures.h"

#include "CAFAna/pca/genie_diag_utils.h"
#include "CAFAna/pca/genie_plot_utils.h"

#include "TCanvas.h"
#include "TH1.h"
#include "TH1D.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 "TLatex.h"
#include "TGaxis.h"

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

using namespace ana;

const int nuniverses = 1000;

// store flavor spectra and their universes
std::vector<TH1D*> nom_hist;
std::vector<TH1D*> shifts_plus;
std::vector<TH1D*> shifts_minus;
std::vector<std::vector<TH1D*>> univ_hist;
  
// 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;

void FillHists(std::string fileName, std::vector<TString> samples) 
{
  for(unsigned int i = 0;i < samples.size();++i){
    // This doesn't actually matter since we'll be normalizing
    double pot = samples[i].Contains("FHC") ? kAna2018FHCPOT : kAna2018RHCPOT;

    // Load the spectra
    std::vector<TH1D*> hUni;
    for(int j = 0;j < nuniverses;++j)
      hUni.push_back(LoadFromFile<Spectrum>(fileName,("genieverse_"+(std::string)samples[i]+"/Univ"+std::to_string(j)).c_str()).release()->ToTH1(pot));

    nom_hist.push_back(hUni[0]);
    
    shifts_plus.push_back(ErrorBand(hUni[0],hUni,1));
    shifts_minus.push_back(ErrorBand(hUni[0],hUni,-1));

    std::vector<TH1D*> temp_hist;
    for(unsigned int j = 1;j < hUni.size();++j)
      temp_hist.push_back(hUni[j]);
    univ_hist.push_back(temp_hist);
  }

}

void goFNBasis(){
  unsigned int length = nom_hist.size();

  // Divide F/N for each component
  for(unsigned int i = length/2;i < length;++i){
    // Numu doesn't extrapolate NCs
    if(i == length/2+4 || i == length/2+5) continue;
    // Nue signal in the ND is numus
    if(i == length/2+6 || i == length/2+7 || i == length/2+8 || i == length/2+9){
      nom_hist[i]->Divide(nom_hist[i-length/2-6]);
      shifts_plus[i]->Divide(shifts_plus[i-length/2-6]);
      shifts_minus[i]->Divide(shifts_minus[i-length/2-6]);
      continue;
    }
    // Others are standard ratios
    nom_hist[i]->Divide(nom_hist[i-length/2]);
    shifts_plus[i]->Divide(shifts_plus[i-length/2]);
    shifts_minus[i]->Divide(shifts_minus[i-length/2]);
  }

  // Same in each universe
  for(unsigned int j = 0;j < univ_hist[0].size();++j){
    for(unsigned int i = length/2;i < length;++i){
      if(i == length/2+4 || i == length/2+5) continue;
      if(i == length/2+6 || i == length/2+7 || i == length/2+8 || i == length/2+9){
        univ_hist[i][j]->Divide(univ_hist[i-length/2-6][j]);
        continue;
      }
      univ_hist[i][j]->Divide(univ_hist[i-length/2][j]);
    }
  }

}

void SuppressND(double sup){
  unsigned int length = nom_hist.size();

  for(unsigned int i = 0;i < length/2;++i){
    nom_hist[i]->Scale(1/sup);
    shifts_plus[i]->Scale(1/sup);
    shifts_minus[i]->Scale(1/sup);
    
    for(unsigned int j = 0;j < univ_hist[0].size();++j)
      univ_hist[i][j]->Scale(1/sup);
  }

}

void NormalizeHists(){
  for(unsigned int i = 0;i < nom_hist.size();++i){
    double nom_area = nom_hist[i]->Integral();
    if(nom_area == 0) continue;

    nom_hist[i]->Scale(1/nom_area);
    shifts_plus[i]->Scale(1/nom_area);
    shifts_minus[i]->Scale(1/nom_area);

    for(unsigned int j = 0;j < univ_hist[i].size();++j)
      univ_hist[i][j]->Scale(1/nom_area);
  }
}

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

    TVectorD x = eigenvectors[component];

    x *= sqrt(eigenvalues[component]);

    TH1D* retplus = new TH1D(x);
    std::string nameplus = "pcplus"+std::to_string(component);
    retplus->SetName(nameplus.c_str());
    retplus->Add(joint_hist_nom);
    pc_shifts_plus.push_back(retplus);

    x *= -1.;

    TH1D* retminus = new TH1D(x);
    std::string nameminus = "pcminus"+std::to_string(component);
    retminus->SetName(nameminus.c_str());
    retminus->Add(joint_hist_nom);
    pc_shifts_minus.push_back(retminus);
    
    std::vector<TH1D*> h_proj_univs;
    // These take /forever/ to calculate with a lot of universes
    for(unsigned int i = 0; i < univ_hist[0].size(); i+=10){

      double projection = 0;

      for(int binIdx = 1; binIdx <= joint_hist_systs[i]->GetNbinsX(); binIdx++) 
        projection += -1.*x[binIdx-1]*(joint_hist_systs[i]->GetBinContent(binIdx)-joint_hist_nom->GetBinContent(binIdx));
     
      projection *= 1./(x.Norm2Sqr());

      TVectorD y = x;
      y *= projection;

      TH1D* ret = new TH1D(y);
      std::string name = "proj"+std::to_string(component)+"univ"+std::to_string(i);
      ret->SetName(name.c_str());
      ret->Add(joint_hist_nom);
      h_proj_univs.push_back(ret);

    }
    pc_projections.push_back(h_proj_univs);
    
  }

}

void SavePCAShifts(int ncomponents,
                   TH1D* joint_hist_nom, std::vector<TH1D*> joint_hist_systs)
{
  TFile* saveHists = new TFile("pc_shifts.root", "RECREATE");
  TDirectory *nomdir = saveHists->mkdir("nom");
  nomdir->cd();
  TH1D* nom_save_hist = (TH1D*)joint_hist_nom->Clone("genie_nom");
  nom_save_hist->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 = "genie_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* genieplus = (TH1D*)pc_shift->Clone("");
    std::string nameplus_save = "genie"+std::to_string(pcIdx);
    genieplus->SetName(nameplus_save.c_str());
    genieplus->Write();
    if(pcIdx >= ncomponents) break;
    pcIdx++;
  }
  saveHists->cd();
  saveHists->Close();
}

void genie_syst_pca()
{

  std::vector<TString> samples = {
    "FHC_NumuND_numus_QE",
    "RHC_NumuND_numus_QE",
    "FHC_NumuND_numus_nonQE",
    "RHC_NumuND_numus_nonQE",
    "FHC_NumuND_ncs_All",
    "RHC_NumuND_ncs_All",

    "FHC_NueND_nues_QE",
    "RHC_NueND_nues_QE",
    "FHC_NueND_nues_nonQE",
    "RHC_NueND_nues_nonQE",
    "FHC_NueND_numus_All",
    "RHC_NueND_numus_All",
    "FHC_NueND_ncs_All",
    "RHC_NueND_ncs_All",

    "FHC_NumuFD_numus_QE",
    "RHC_NumuFD_numus_QE",
    "FHC_NumuFD_numus_nonQE",
    "RHC_NumuFD_numus_nonQE",
    "FHC_NumuFD_ncs_All",
    "RHC_NumuFD_ncs_All",

    "FHC_NueFD_nues_QE",
    "RHC_NueFD_nues_QE",
    "FHC_NueFD_nues_nonQE",
    "RHC_NueFD_nues_nonQE",
    "FHC_NueFD_numus_All",
    "RHC_NueFD_numus_All",
    "FHC_NueFD_ncs_All",
    "RHC_NueFD_ncs_All",
  };  

  // input root file      
  std::string fileName = "make_genieverse.root";
  
  // load the spectra
  FillHists(fileName, samples);
  // Switch to F/N for the relevant spectra
  goFNBasis();
  // Normalize importance
  NormalizeHists();
  // Suppress ND to improve convergence
  SuppressND(suppress);
  
  // flatten spectra
  TH1D* joint_hist_nom = CombineHistograms(nom_hist);
  TH1D* joint_hist_plus = CombineHistograms(shifts_plus);
  TH1D* joint_hist_minus = CombineHistograms(shifts_minus);
  std::vector <TH1D*> joint_hist_systs;
  for(unsigned int i = 0;i < univ_hist[0].size();++i){
    std::vector<TH1D*> temp_hist;
    for(unsigned int j = 0;j < univ_hist.size();++j) temp_hist.push_back(univ_hist[j][i]);
    TH1D* joint_hist_univ = CombineHistograms(temp_hist);
    joint_hist_systs.push_back(joint_hist_univ);
  }

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

  // calculate bin-bin covariance
  const std::vector<TH1D*> genie(joint_hist_systs.begin(), joint_hist_systs.end());
  std::unique_ptr<TMatrixDSym> covMx = GetCovMx(genie);
  
  // 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);

  // Get PCA Shifts and Projections
  FillPCAContainers(eigenvectors, eigenvalues, joint_hist_nom, joint_hist_systs);  

  // Save Shifts
  int ncomponents = 100;
  SavePCAShifts(ncomponents, joint_hist_nom, joint_hist_systs);

  // plotting!  
  // plot the distribution of universes
  joint_hist_nom->SetMaximum(1.4*joint_hist_nom->GetMaximum());
  joint_hist_nom->GetXaxis()->CenterTitle();
  joint_hist_plus->SetLineColor(kAzure);
  joint_hist_minus->SetLineColor(kRed);
  std::vector<TH1D*> joints;
  for(unsigned int i = 0;i < joint_hist_systs.size();++i){
    joint_hist_systs[i]->SetLineColor(kGray);
    joints.push_back(joint_hist_systs[i]);
  }
  joints.push_back(joint_hist_minus);
  joints.push_back(joint_hist_plus);

  TCanvas *c1 = QuickUnivRatioPlot("c1",joint_hist_nom,joints);

  TLegend *l = new TLegend(0.64,0.69,0.92,0.84);
  l->SetFillStyle(0);
  l->AddEntry(joint_hist_nom,"Nominal MC","l");
  l->AddEntry(joint_hist_plus,"+1#sigma Universe","l");
  l->AddEntry(joint_hist_minus,"-1#sigma Universe","l");
  l->AddEntry(joint_hist_systs[0],"Genie Universe","l");
  l->Draw();

  drawLabels(axis_labels);

  c1->cd(1);
  drawBinLines(axis_widths);
  c1->cd(2);
  drawBinLines(axis_widths);

  c1->Print("plots/genie_all.pdf");

  // plot the covariance matrix
  TH2D *histCov = new TH2D("histCov", "Bin-Bin Covariance", joint_nbins, 0, joint_nbins, joint_nbins, 0, joint_nbins);
  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]);
    }
  }

  gStyle->SetPalette(kCool);

  TCanvas *c2 = new TCanvas("c2", "c2");
  histCov->GetXaxis()->CenterTitle();
  histCov->GetYaxis()->CenterTitle();
  histCov->Draw("COLZ");
  c2->SetLogz();
  c2->SetRightMargin(0.125);
  c2->Update();

  //drawLabels2D(axis_labels);

  c2->Print("plots/genie_covariance.pdf");

  for(int i = 0;i < eigenvalues.GetNrows();++i)
    if(eigenvalues[i] < 10E-15){
      std::cout<<i-2<<" "<<eigenvalues[i-2]<<std::endl;
      std::cout<<i-1<<" "<<eigenvalues[i-1]<<std::endl;
      std::cout<<i  <<" "<<eigenvalues[i]  <<std::endl;
      std::cout<<i+1<<" "<<eigenvalues[i+1]<<std::endl;
      std::cout<<i+2<<" "<<eigenvalues[i+2]<<std::endl;
      break;
    }

  // plot the eigenvalue profile
  TCanvas *c3 = new TCanvas("c3", "c3");
  TH1D* histEigen = new TH1D(eigenvalues);
  histEigen->SetName("pi");
  histEigen->GetXaxis()->SetTitle("Eigennumber");
  histEigen->GetYaxis()->SetTitle("Eigenvalue");
  histEigen->SetMarkerStyle(kFullCircle);
  histEigen->GetXaxis()->SetRangeUser(0,238);
  histEigen->GetXaxis()->CenterTitle();
  histEigen->GetYaxis()->CenterTitle();
  histEigen->Draw("p");
  c3->SetLogy();

  c3->Print("plots/genie_eigen.pdf");

  // plot the total explained variancee
  TCanvas *c4 = new TCanvas("c4","c4");
  TVectorD eigenInt(eigenvalues.GetNrows());
  double sum=0;
  for(int i = 0;i < eigenvalues.GetNrows();++i){
    sum+=eigenvalues[i];
    eigenInt[i] = sum/eigenvalues.Sum();
  }
  TH1D* histEigenInt = new TH1D(eigenInt);
  histEigenInt->GetXaxis()->SetTitle("N_PC");
  histEigenInt->GetYaxis()->SetTitle("Total Explained Variance");
  histEigenInt->SetMarkerStyle(kFullCircle);
  histEigenInt->GetXaxis()->SetRangeUser(0,238);
  histEigenInt->GetXaxis()->CenterTitle();
  histEigenInt->GetYaxis()->CenterTitle();
  histEigenInt->Draw("p");

  c4->Update();
  drawCoverLines(histEigenInt);

  c4->Print("plots/genie_coverage.pdf");
  
  // plot the PC and the universe projections
  //for(int i = 0;i < eigenvalues.GetNrows();++i){
  // only like 20 tho... No one will ever look at more than that right?
  for(int i = 0;i < 20;++i){
    for(unsigned int j = 0;j < pc_projections[i].size();++j) pc_projections[i][j]->SetLineColor(kGray);
    pc_shifts_plus[i]->SetLineColor(kAzure);
    pc_shifts_minus[i]->SetLineColor(kRed);
    std::vector<TH1D*> hComb(pc_projections[i].begin(),pc_projections[i].end());
    hComb.push_back(pc_shifts_minus[i]);
    hComb.push_back(pc_shifts_plus[i]);
    TCanvas *cproj = QuickUnivRatioPlot("cproj_"+std::to_string(i),joint_hist_nom,hComb);    

    drawLabels(axis_labels);

    TLegend *lproj = new TLegend(0.64,0.69,0.92,0.84);
    lproj->SetFillStyle(0);
    lproj->AddEntry(joint_hist_nom,"Nominal MC","l");
    lproj->AddEntry(pc_shifts_plus[i],"+1#sigma Shift","l");
    lproj->AddEntry(pc_shifts_minus[i],"-1#sigma Shift","l");
    lproj->AddEntry(pc_projections[i][0],"Universes Projections","l");
    lproj->Draw();

    cproj->cd(1);
    drawBinLines(axis_widths);
    cproj->cd(2);
    drawBinLines(axis_widths);

    cproj->cd();
    TLatex *tex = new TLatex(0.9,0.93,("PC"+std::to_string(i)+" Explained Variance: "+std::to_string(eigenvalues[i]/eigenvalues.Sum())).c_str());
    tex->SetTextSize(1./30);
    tex->SetTextAlign(32);
    tex->Draw();
    
    if(i<10) cproj->Print(("plots/genie_proj_PC00"+std::to_string(i)+".pdf").c_str());
    else if(i<100) cproj->Print(("plots/genie_proj_PC0"+std::to_string(i)+".pdf").c_str());
    else std::cout<<"Uh, what?"<<std::endl;
    
  }

}