Functions | Variables
ppfx_syst_pca_fn.C File Reference
#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>

Go to the source code of this file.

Functions

void FillFlavorContainers (std::string fileName, std::vector< std::string > comps, std::vector< std::string > comps2)
 
TH1D * ConvertNormalBasis (TH1D *joint_hist_fn)
 
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)
 
void SavePCAShifts (int ncomponents, TH1D *joint_hist_cv, std::vector< TH1D * > joint_hist_systs)
 
void SavePCAForCAFAna (int ncomponents, double inflate=1.)
 
TH1D * GetFNRatio (TH1D *joint_hist)
 
TH1D * GetTotFNRatio (TH1D *joint_hist)
 
TH1D * FNBasisHelper (TH1D *near, TH1D *far, double near_norm, double far_norm)
 
void ConvertFNBasis (double nonfn_scale=1e-6)
 
void ppfx_syst_pca_fn ()
 

Variables

std::vector< std::string > comps
 
std::vector< std::string > comps2
 
int rebin_f = 10
 
int nunivs = 100
 
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
 
std::vector< TH1D * > pc_shifts_plus
 
std::vector< TH1D * > pc_shifts_minus
 
TH1D * fn_cv
 
std::vector< TH1D * > fn_shifts_ppfx
 
std::vector< TH1D * > fn_shifts_plus
 
std::vector< TH1D * > fn_shifts_minus
 

Function Documentation

void ConvertFNBasis ( double  nonfn_scale = 1e-6)

Definition at line 352 of file ppfx_syst_pca_fn.C.

References plot_validation_datamc::Clone(), cv_areas, cv_areas_fn, cv_comps, cv_comps_fn, FNBasisHelper(), makeTrainCVSamples::int, nunivs, ana::UniqueName(), univs_comps, and univs_comps_fn.

Referenced by ppfx_syst_pca_fn().

352  {
353 
354  for(int nearIdx = 0; nearIdx < 8; nearIdx++) {
355  TH1D* near = (TH1D*)cv_comps[nearIdx]->Clone(UniqueName().c_str());
356  near->Scale(nonfn_scale);
357  cv_comps_fn.push_back(near);
358  cv_areas_fn.push_back(cv_areas[nearIdx]/nonfn_scale);
359  }
360  for(int farIdx = 8; farIdx < 16; farIdx++) {
361  TH1D* fn = FNBasisHelper(cv_comps[farIdx-8], cv_comps[farIdx], cv_areas[farIdx-8], cv_areas[farIdx]);
362  cv_areas_fn.push_back(fn->Integral());
363  fn->Scale(1./fn->Integral());
364  cv_comps_fn.push_back(fn);
365  }
366 
367  for(int UnivIdx = 0; UnivIdx < (int)nunivs; UnivIdx++){
368  std::vector<TH1D*> temp_fn_univ;
369  for(int nearIdx = 0; nearIdx < 8; nearIdx++) {
370  TH1D* near_univ = (TH1D*)univs_comps[UnivIdx][nearIdx]->Clone(UniqueName().c_str());
371  near_univ->Scale(nonfn_scale);
372  temp_fn_univ.push_back(near_univ);
373  }
374  for(int farIdx = 8; farIdx < 16; farIdx++) {
375  TH1D* fn_univ = FNBasisHelper(univs_comps[UnivIdx][farIdx-8], univs_comps[UnivIdx][farIdx], cv_areas[farIdx-8], cv_areas[farIdx]);
376  fn_univ->Scale(1./cv_areas_fn[farIdx]);
377  temp_fn_univ.push_back(fn_univ);
378  }
379  univs_comps_fn.push_back(temp_fn_univ);
380  }
381 }
TH1D * FNBasisHelper(TH1D *near, TH1D *far, double near_norm, double far_norm)
std::vector< double > cv_areas
std::vector< TH1D * > cv_comps_fn
std::vector< TH1D * > cv_comps
std::vector< std::vector< TH1D * > > univs_comps
std::vector< std::vector< TH1D * > > univs_comps_fn
std::vector< double > cv_areas_fn
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:30
int nunivs
TH1D* ConvertNormalBasis ( TH1D *  joint_hist_fn)

Definition at line 121 of file ppfx_syst_pca_fn.C.

References plot_validation_datamc::Clone(), CombineHistograms(), comps, cv_areas, cv_areas_fn, makeTrainCVSamples::int, runNovaSAM::ret, SplitHistograms(), and ana::UniqueName().

Referenced by FillPCAContainers(), GetTotFNRatio(), and ppfx_syst_pca_fn().

121  {
122 
123  std::vector<TH1D*> fn_basis = SplitHistograms(joint_hist_fn, comps.size());
124  std::vector<TH1D*> near_hists(fn_basis.begin(), fn_basis.begin()+8);
125  std::vector<TH1D*> fn_hists(fn_basis.begin()+8, fn_basis.end());
126 
127  std::vector<TH1D*> normal_basis;
128  for(int nearIdx = 0; nearIdx < (int)near_hists.size(); nearIdx++){
129  TH1D* h_normal = (TH1D*)near_hists[nearIdx]->Clone(UniqueName().c_str());
130  h_normal->Scale(cv_areas_fn[nearIdx]/cv_areas[nearIdx]);
131  normal_basis.push_back(h_normal);
132  }
133  for(int fnIdx = 0; fnIdx < (int)fn_hists.size(); fnIdx++){
134  TH1D* far_h = (TH1D*)fn_hists[fnIdx]->Clone(UniqueName().c_str());
135  far_h->Scale(cv_areas_fn[fnIdx+8]);
136  far_h->Multiply(normal_basis[fnIdx]);
137  far_h->Scale(cv_areas[fnIdx]/cv_areas[fnIdx+8]);
138  normal_basis.push_back(far_h);
139  }
140  TH1D* ret = CombineHistograms(normal_basis);
141  return ret;
142 }
std::vector< TH1D * > SplitHistograms(TH1D *hist, int nhists)
TH1D * CombineHistograms(std::vector< TH1D * > hists)
std::vector< double > cv_areas
std::vector< std::string > comps
std::vector< double > cv_areas_fn
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:30
void FillFlavorContainers ( std::string  fileName,
std::vector< std::string >  comps,
std::vector< std::string >  comps2 
)

Definition at line 88 of file ppfx_syst_pca_fn.C.

References cv_areas, cv_comps, genie::utils::style::Format(), makeTrainCVSamples::int, ana::Spectrum::LoadFrom(), nunivs, ana::rebin_f, and univs_comps.

Referenced by ppfx_syst_pca_fn().

89 {
90  TFile* fload = new TFile(fileName.c_str(), "READ");
91 
92  for(int compIdx = 0; compIdx < (int)comps.size(); ++compIdx){
93  const char* comp = comps[compIdx].c_str();
94 
95  auto spec = Spectrum::LoadFrom(fload, TString::Format("spect_%s", comp));
96 
97  TH1D* temp_hist = (TH1D*)spec->ToTH1(spec->POT());
98  temp_hist->Rebin(rebin_f);
99  cv_areas.push_back(temp_hist->Integral());
100  temp_hist->Scale(1./temp_hist->Integral());
101  cv_comps.push_back(temp_hist);
102  }
103 
104  // universes are normalized to cv areas
105  for(int UnivIdx = 0; UnivIdx < nunivs; ++UnivIdx){
106  std::vector<TH1D*> temp_hist_systs;
107  for(int compIdx = 0; compIdx < (int)comps2.size(); ++compIdx){
108  const char* comp = comps2[compIdx].c_str();
109  auto spec_syst = Spectrum::LoadFrom(fload, TString::Format("spect_%s_%d", comp, UnivIdx));
110 
111  TH1D* hist_syst = (TH1D*)spec_syst->ToTH1(spec_syst->POT());
112  hist_syst->Rebin(rebin_f);
113  hist_syst->Scale(1./cv_areas[compIdx]);
114  temp_hist_systs.push_back(hist_syst);
115 
116  }
117  univs_comps.push_back(temp_hist_systs);
118  }
119 }
std::vector< std::string > comps2
fileName
Definition: plotROC.py:78
int rebin_f
std::vector< double > cv_areas
std::unique_ptr< T > LoadFrom(TDirectory *dir, const std::string &label)
Definition: LoadFromFile.h:17
std::vector< TH1D * > cv_comps
std::vector< std::vector< TH1D * > > univs_comps
std::vector< std::string > comps
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
int nunivs
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 
)

Definition at line 144 of file ppfx_syst_pca_fn.C.

References test_component::component, ConvertNormalBasis(), pc_shifts_minus, pc_shifts_plus, std::sqrt(), art::to_string(), ana::UniqueName(), and submit_syst::x.

Referenced by ppfx_syst_pca_fn().

147 {
148  for(int component = 0; component < eigenvalues.GetNrows(); component++){
149 
150  TVectorD x = eigenvectors[component];
151 
152  x *= sqrt(eigenvalues[component]);
153 
154  TH1D* retplus_fn = new TH1D(x);
155  retplus_fn->SetName(UniqueName().c_str());
156  retplus_fn->Add(joint_hist_cv_fn);
157  TH1D* retplus = ConvertNormalBasis(retplus_fn);
158  std::string nameplus = "pcplus"+std::to_string(component);
159  retplus->SetName(nameplus.c_str());
160  retplus->Divide(joint_hist_cv);
161  retplus->GetYaxis()->SetRangeUser(0.7, 1.3);
162  pc_shifts_plus.push_back(retplus);
163 
164  x *= -1.;
165 
166  TH1D* retminus_fn = new TH1D(x);
167  retminus_fn->SetName(UniqueName().c_str());
168  retminus_fn->Add(joint_hist_cv_fn);
169  TH1D* retminus = ConvertNormalBasis(retminus_fn);
170  std::string nameminus = "pcminus"+std::to_string(component);
171  retminus->SetName(nameminus.c_str());
172  retminus->Divide(joint_hist_cv); // we only care about fractions
173  retminus->GetYaxis()->SetRangeUser(0.7, 1.3);
174  pc_shifts_minus.push_back(retminus);
175 
176  }
177 
178 }
T sqrt(T number)
Definition: d0nt_math.hpp:156
std::vector< TH1D * > pc_shifts_minus
std::vector< TH1D * > pc_shifts_plus
TH1D * ConvertNormalBasis(TH1D *joint_hist_fn)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:30
TH1D* FNBasisHelper ( TH1D *  near,
TH1D *  far,
double  near_norm,
double  far_norm 
)

Definition at line 343 of file ppfx_syst_pca_fn.C.

References ana::UniqueName().

Referenced by ConvertFNBasis().

344 {
345  TH1D* fn = (TH1D*)far->Clone(UniqueName().c_str());
346  fn->Divide(near);
347  fn->Scale(far_norm/near_norm);
348 
349  return fn;
350 }
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:30
TH1D* GetFNRatio ( TH1D *  joint_hist)

Definition at line 312 of file ppfx_syst_pca_fn.C.

References plot_validation_datamc::Clone(), CombineHistograms(), comps, cv_areas, makeTrainCVSamples::int, runNovaSAM::ret, Scale(), SplitHistograms(), and ana::UniqueName().

Referenced by GetTotFNRatio(), and ppfx_syst_pca_fn().

313 {
314  int ncomps = (int) comps.size();
315  std::vector<TH1D*> hflavs = SplitHistograms(joint_hist, ncomps);
316  for(int comp = 0; comp < ncomps; comp++)
317  hflavs[comp]->Scale(cv_areas[comp]);
318 
319  std::vector<TH1D*> near_flavs(hflavs.begin(), hflavs.begin()+8);
320  std::vector<TH1D*> far_flavs(hflavs.begin()+8, hflavs.end());
321  std::vector<TH1D*> fn_flavs;
322  for(int flavIdx = 0; flavIdx < (int) near_flavs.size(); flavIdx++){
323  TH1D* hfn = (TH1D*)far_flavs[flavIdx]->Clone(UniqueName().c_str());
324  hfn->Divide(near_flavs[flavIdx]);
325  fn_flavs.push_back(hfn);
326  }
327  TH1D* ret = CombineHistograms(fn_flavs);
328 
329  return ret;
330 
331 }
std::vector< TH1D * > SplitHistograms(TH1D *hist, int nhists)
TH1D * CombineHistograms(std::vector< TH1D * > hists)
std::vector< double > cv_areas
std::vector< std::string > comps
simulatedPeak Scale(1/simulationNormalisationFactor)
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:30
TH1D* GetTotFNRatio ( TH1D *  joint_hist)

Definition at line 333 of file ppfx_syst_pca_fn.C.

References ConvertNormalBasis(), GetFNRatio(), analysePickle::hist, and runNovaSAM::ret.

Referenced by ppfx_syst_pca_fn().

334 {
335 
336  TH1D* hist = ConvertNormalBasis(joint_hist);
337  TH1D* ret = GetFNRatio(hist);
338 
339  return ret;
340 
341 }
TH1D * GetFNRatio(TH1D *joint_hist)
TH1D * ConvertNormalBasis(TH1D *joint_hist_fn)
void ppfx_syst_pca_fn ( )

Definition at line 384 of file ppfx_syst_pca_fn.C.

References axis_labels, demo5::c1, demo5::c2, chisquared::c3, chisquared::c4, plot_validation_datamc::Clone(), CombineHistograms(), test_component::component, comps, comps2, ConvertFNBasis(), ConvertNormalBasis(), om::cout, CrossCheckDiag(), cv, cv_comps_fn, draw_vlines(), plotROC::fileName, FillFlavorContainers(), FillPCAContainers(), fn_cv, fn_shifts_minus, fn_shifts_plus, fn_shifts_ppfx, GetCovMx(), GetFNRatio(), GetTotFNRatio(), make_syst_table_plots::h, analysePickle::hist, MECModelEnuComparisons::i, make_pair(), nbins, nunivs, pc_shifts_minus, pc_shifts_plus, plot_hists(), plot_label(), plot_ratios(), SavePCAForCAFAna(), SavePCAShifts(), set_labels_1D(), set_labels_2D(), art::to_string(), ana::UniqueName(), and univs_comps_fn.

385 {
386 
387  // input root file
388  std::string fileName = "ppfx_smooth_flux.root";
389 
390  // load the spectra and normalize
391  FillFlavorContainers(fileName, comps, comps2);
392  ConvertFNBasis();
393 
394  // flatten spectra
395  TH1D* joint_hist_cv_fn = CombineHistograms(cv_comps_fn);
396  std::vector <TH1D*> joint_hist_systs_fn;
397  for(auto hist: univs_comps_fn) {
398  TH1D* joint_hist_univ = CombineHistograms(hist);
399  joint_hist_systs_fn.push_back(joint_hist_univ);
400  }
401 
402  // get number of bins
403  int nbins = cv_comps_fn[0]->GetNbinsX();
404  int joint_nbins = joint_hist_cv_fn->GetNbinsX();
405 
406  // calculate bin-bin covariance
407  const std::vector<TH1D*> ppfx(joint_hist_systs_fn.begin(), joint_hist_systs_fn.end());
408  std::cout << "nunivs: " << ppfx.size();
409  // std::unique_ptr<TMatrixDSym> covMx = GetCovMx(ppfx);
410  std::unique_ptr<TMatrixDSym> covMx = GetCovMx(ppfx, joint_hist_cv_fn);
411 
412  // diagonalize
413  TMatrixDSymEigen cov(*covMx.get());
414  TMatrixD V = cov.GetEigenVectors();
415  TMatrixD eigenvectors = V;
416  eigenvectors.Transpose(V);
417  TVectorD eigenvalues = cov.GetEigenValues();
418 
419  // Verify that the diagonalisation works as expected
420  CrossCheckDiag(*covMx.get(), V, eigenvectors, eigenvalues);
421 
422  TH1D* joint_hist_cv = ConvertNormalBasis(joint_hist_cv_fn);
423  std::vector<TH1D*> joint_hist_systs;
424  std::for_each(joint_hist_systs_fn.begin(), joint_hist_systs_fn.end(), [&](TH1D* joint_h){
425  joint_hist_systs.push_back(ConvertNormalBasis(joint_h));
426  });
427 
428  // Get PCA Shifts and Projections
429  FillPCAContainers(eigenvectors, eigenvalues, joint_hist_cv_fn, joint_hist_systs_fn, joint_hist_cv, joint_hist_systs);
430 
431  // store only these principle components
432  int ncomponents = 100;
433 
434  // sum components up to get far/near ratios
435  fn_cv = GetTotFNRatio(joint_hist_cv_fn);
436  for(int i = 0; i < nunivs; i++)
437  fn_shifts_ppfx.push_back(GetTotFNRatio(joint_hist_systs_fn[i]));
438 
439  for(int i = 0; i < ncomponents; i++){
440  TH1D* fn_plus = (TH1D*)pc_shifts_plus[i]->Clone(UniqueName().c_str());
441  fn_plus->Multiply(joint_hist_cv);
442  fn_shifts_plus.push_back(GetFNRatio(fn_plus));
443 
444  TH1D* fn_minus = (TH1D*)pc_shifts_minus[i]->Clone(UniqueName().c_str());
445  fn_minus->Multiply(joint_hist_cv);
446  fn_shifts_minus.push_back(GetFNRatio(fn_minus));
447  }
448 
449  SavePCAShifts(ncomponents, joint_hist_cv, joint_hist_systs);
450  SavePCAForCAFAna(50, 1.25);
451 
452  // plotting!
453  // store axis labels
454 
455  std::vector<std::string> axis_labels = {
456  "#nu_{#mu} ND RHC", "#nu_{e} ND RHC", "#bar{#nu}_{#mu} ND RHC", "#bar{#nu}_{e} ND RHC",
457  "#nu_{#mu} ND FHC", "#nu_{e} ND FHC", "#bar{#nu}_{#mu} ND FHC", "#bar{#nu}_{e} ND FHC",
458  "#nu_{#mu} F/N RHC", "#nu_{e} F/N RHC", "#bar{#nu}_{#mu} F/N RHC", "#bar{#nu}_{e} F/N RHC",
459  "#nu_{#mu} F/N FHC", "#nu_{e} F/N FHC", "#bar{#nu}_{#mu} F/N FHC", "#bar{#nu}_{e} F/N FHC",
460  };
461  std::vector<std::string> axis_labels_fn = {
462  "#nu_{#mu} F/N RHC", "#nu_{e} F/N RHC", "#bar{#nu}_{#mu} F/N RHC", "#bar{#nu}_{e} F/N RHC",
463  "#nu_{#mu} F/N FHC", "#nu_{e} F/N FHC", "#bar{#nu}_{#mu} F/N FHC", "#bar{#nu}_{e} F/N FHC",
464  };
465  std::vector<std::string> axis_labels_nom = {
466  "#nu_{#mu} ND RHC", "#nu_{e} ND RHC", "#bar{#nu}_{#mu} ND RHC", "#bar{#nu}_{e} ND RHC",
467  "#nu_{#mu} ND FHC", "#nu_{e} ND FHC", "#bar{#nu}_{#mu} ND FHC", "#bar{#nu}_{e} ND FHC",
468  "#nu_{#mu} FD RHC", "#nu_{e} FD RHC", "#bar{#nu}_{#mu} FD RHC", "#bar{#nu}_{e} FD RHC",
469  "#nu_{#mu} FD FHC", "#nu_{e} FD FHC", "#bar{#nu}_{#mu} FD FHC", "#bar{#nu}_{e} FD FHC",
470  };
471 
472  // custom methods for pre-diagonalization plotting
473  std::vector <plot_h> h_cv_systs;
474  std::vector <std::pair<plot_h, TH1*>> h_cv_systs_ratios;
475 
476  for(auto hist: joint_hist_systs){
477  plot_h univ = {hist, "hist same", kRed - 7, false, "", "l"};
478  h_cv_systs.push_back(univ);
479  h_cv_systs_ratios.push_back(std::make_pair(univ, joint_hist_cv));
480  }
481 
482  plot_h cv = {joint_hist_cv, "hist same", kBlack, false, "", "l"};
483  h_cv_systs.push_back(cv);
484  h_cv_systs_ratios.push_back(std::make_pair(cv, joint_hist_cv));
485 
486  TCanvas *c1 = new TCanvas("c1", "c1", 200, 10, 800, 600);
487  c1->cd();
488  c1->Draw();
489  plot_hists(h_cv_systs, "", false, false);
490  set_labels_1D(axis_labels_nom, joint_hist_cv);
491  draw_vlines(comps.size(), joint_hist_cv->GetXaxis()->GetBinLowEdge(joint_nbins+1), 0, 0.26, kBlue);
492  c1->Update();
493  c1->Print("ppfx_diag_plots/ppfx_hadp_spectra.pdf");
494 
495  TCanvas *c2 = new TCanvas("c2", "c2", 200, 10, 800, 600);
496  c2->cd();
497  c2->Draw();
498  plot_ratios(h_cv_systs_ratios, "", false, 0, 2);
499  set_labels_1D(axis_labels_nom, joint_hist_cv);
500  draw_vlines(comps.size(), joint_hist_cv->GetXaxis()->GetBinLowEdge(joint_nbins+1), 0, 2);
501  c2->Update();
502  c2->Print("ppfx_diag_plots/ppfx_hadp_spectra_ratios.pdf");
503 
504  // plot the covariance matrix
505  TH2D *histCov = new TH2D("histCov", "Covariance", joint_nbins, 0, joint_nbins, joint_nbins, 0, joint_nbins);
506  TH2D *histCov_fn = new TH2D("histCov_fn", "Covariance", joint_nbins/2, 0, joint_nbins/2, joint_nbins/2, 0, joint_nbins/2);
507  for(int xbinIdx = 1; xbinIdx <= joint_nbins; xbinIdx++){
508  for(int ybinIdx = 1; ybinIdx <= joint_nbins; ybinIdx++){
509  histCov->Fill(xbinIdx, ybinIdx, (*covMx)[xbinIdx-1][ybinIdx-1]);
510  }
511  }
512  for(int xbinIdx = (joint_nbins/2)+1; xbinIdx <= joint_nbins; xbinIdx++){
513  for(int ybinIdx = (joint_nbins/2)+1; ybinIdx <= joint_nbins; ybinIdx++){
514  histCov_fn->Fill(xbinIdx-(joint_nbins/2), ybinIdx-(joint_nbins/2), (*covMx)[xbinIdx-1][ybinIdx-1]);
515  }
516  }
517 
518  gStyle->SetPalette(55);
519  TCanvas *c3 = new TCanvas("c3", "c3", 200, 10, 800, 600);
520  c3->cd();
521  c3->Draw();
522  histCov->Draw("COLZ");
523  c3->SetLogz();
524  c3->SetRightMargin(0.125);
525  set_labels_2D(axis_labels, axis_labels, histCov);
526  c3->Update();
527  c3->Print("ppfx_diag_plots/ppfx_hadp_covariance.pdf");
528 
529  TCanvas *c3_fn = new TCanvas("c3_fn", "c3_fn", 200, 10, 800, 600);
530  c3_fn->cd();
531  c3_fn->Draw();
532  histCov_fn->Draw("COLZ");
533  c3_fn->SetLogz();
534  c3_fn->SetRightMargin(0.125);
535  set_labels_2D(axis_labels_fn, axis_labels_fn, histCov_fn);
536  c3_fn->Update();
537  c3_fn->Print("ppfx_diag_plots/ppfx_hadp_covariance_fn.pdf");
538 
539  // plot the eigenvalue profile
540  TH1D* histEigen = new TH1D(eigenvalues);
541  histEigen->GetXaxis()->SetTitle("Eigennumber");
542  histEigen->GetYaxis()->SetTitle("Eigenvalue");
543  TCanvas *c4 = new TCanvas("c4", "c4", 200, 10, 800, 600);
544  c4->cd();
545  c4->Draw();
546  histEigen->Draw("l");
547  c4->SetLogy();
548  c4->Update();
549  c4->Print("ppfx_diag_plots/ppfx_hadp_eigenvalues.pdf");
550 
551  int nppfx_diag_plots = 5;
552  // plot the far/near ratio and shifts from the central value
553  for(int component = 0; component < nppfx_diag_plots; component++){
554 
555  std::vector<plot_h> far_near_shifts = {
556  {fn_cv, "hist same", kBlack, true, "CV", "l"},
557  {fn_shifts_plus[component], "hist same", kBlue + 1, true, "Shift Plus", "l"},
558  {fn_shifts_minus[component], "hist same", kRed, true, "Shift Minus", "l"},
559  };
560 
561  std::vector<std::pair<plot_h, TH1*>> double_ratios;
562  std::for_each(far_near_shifts.begin(), far_near_shifts.end(),
563  [&](plot_h h){ double_ratios.push_back(std::make_pair(h, fn_cv));
564  });
565 
566  std::vector<plot_h> pc_shifts;
567  // for(int i = 0; i < (int)pc_projections[component].size(); i++)
568  // pc_shifts.push_back({pc_projections[component][i], "hist same", kGray, false, "", "l"});
569 
570  pc_shifts.push_back({pc_shifts_plus[component], "hist same", kBlue + 1, false, "", "l"});
571  pc_shifts.push_back({pc_shifts_minus[component], "hist same", kRed, false, "", "l"});
572 
573  // shifts and multiverse projection ppfx_diag_plots
574  TCanvas *c5 = new TCanvas("c5", "c5", 200, 10, 800, 600);
575  c5->cd();
576  c5->Draw();
577  plot_hists(pc_shifts, "", false, false);
578  plot_label(0.5, 0.93, Form("Proportion of Variance : %lg", eigenvalues[component]/eigenvalues.Sum()));
579  set_labels_1D(axis_labels_nom, pc_shifts_minus[component]);
580  draw_vlines(comps.size(), joint_hist_cv->GetXaxis()->GetBinLowEdge(joint_nbins+1), 0.7, 1.3);
581  c5->Update();
582  std::string filename1 = "ppfx_diag_plots/ppfx_hadp_shifts_pc"+std::to_string(component)+".pdf";
583  c5->Print(filename1.c_str());
584 
585  TCanvas *c6 = new TCanvas("c6", "c6", 200, 10, 800, 600);
586  c6->cd();
587  TPad *padUp = new TPad("padUp","padUp",0.0,0.30,1.,1.00);
588  padUp->Draw();
589  padUp->cd();
590  padUp->SetBottomMargin(0.00);
591  padUp->SetBorderSize(0);
592 
593  // far/near ratio ppfx_diag_plots
594  std::string fn_title = "F/N - Principle Component " + std::to_string(component);
595  plot_hists(far_near_shifts, fn_title, false, false);
596  // TLegend* leg = new TLegend(0.5, 0.15, 0.8, 0.35);
597  // std::for_each(far_near_shifts.begin(), far_near_shifts.end(),
598  // [&](plot_h h){ leg->AddEntry(h.hist, h.legendLabel.c_str(), h.legendOption);
599  // });
600  // leg->Draw();
601  padUp->Update();
602 
603  c6->cd();
604  TPad *padDn = new TPad("padDn","padDn",0.0,0.0,1.,0.30);
605  padDn->Draw();
606  padDn->cd();
607  padDn->SetBottomMargin(0.30);
608  padDn->SetTopMargin(0.00);
609  padDn->SetBorderSize(0);
610  TH1* base = plot_ratios(double_ratios, "", false, 0.975, 1.025);
611  base->GetYaxis()->SetTitle("Double Ratio");
612  base->GetXaxis()->SetTitle("");
613  base->GetYaxis()->SetTitleSize(base->GetYaxis()->GetTitleSize()*7./3.);
614  base->GetYaxis()->SetTitleOffset(base->GetYaxis()->GetTitleOffset()*3./7.);
615  base->GetYaxis()->SetLabelSize(base->GetYaxis()->GetLabelSize()*7./3.);
616  base->GetXaxis()->SetTitleSize(base->GetXaxis()->GetTitleSize()*7./3.);
617  base->GetXaxis()->SetTitleOffset(base->GetXaxis()->GetTitleOffset()*6./7.);
618  base->GetXaxis()->SetLabelSize(base->GetXaxis()->GetLabelSize()*7./3.);
619  set_labels_1D(axis_labels_fn, joint_hist_cv);
620  padDn->Update();
621 
622  std::string filename2 = "ppfx_diag_plots/ppfx_hadp_fn_ratio_pc"+std::to_string(component)+".pdf";
623  c6->Print(filename2.c_str());
624 
625  }
626 }
TH1D * fn_cv
std::vector< std::string > comps2
TH1D * GetTotFNRatio(TH1D *joint_hist)
fileName
Definition: plotROC.py:78
std::vector< TH1D * > fn_shifts_ppfx
TH1D * CombineHistograms(std::vector< TH1D * > hists)
void SavePCAForCAFAna(int ncomponents, double inflate=1.)
std::vector< TH1D * > fn_shifts_minus
void SavePCAShifts(int ncomponents, TH1D *joint_hist_cv, std::vector< TH1D * > joint_hist_systs)
std::vector< TH1D * > pc_shifts_minus
TH1D * GetFNRatio(TH1D *joint_hist)
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
std::vector< TH1D * > pc_shifts_plus
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)
std::vector< TH1D * > cv_comps_fn
void set_labels_2D(std::vector< std::string > labels_x, std::vector< std::string > labels_y, TH2D *hist)
const XML_Char int const XML_Char int const XML_Char * base
Definition: expat.h:331
std::vector< TH1D * > fn_shifts_plus
c2
Definition: demo5.py:33
TH1D * ConvertNormalBasis(TH1D *joint_hist_fn)
const int nbins
Definition: cellShifts.C:15
const std::string cv[Ncv]
void plot_hists(std::vector< plot_h > hists, std::string title="", bool turnOnlegend=true, bool turnOnMaxRange=true, double leg_w=0.2, double leg_h=0.3, double offset=1.5)
std::vector< std::string > comps
OStream cout
Definition: OStream.cxx:6
void ConvertFNBasis(double nonfn_scale=1e-6)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
std::vector< std::vector< TH1D * > > univs_comps_fn
std::unique_ptr< TMatrixDSym > GetCovMx(const std::vector< TH1D * > &hists)
TH1 * plot_ratios(std::vector< std::pair< plot_h, TH1 * >> hists, std::string title="", bool turnOnlegend=false, double low_y=-1, double high_y=-1)
void set_labels_1D(std::vector< std::string > labels, TH1D *hist)
void draw_vlines(int n, double xmax, double y0, double y1, Color_t color=kBlack)
c1
Definition: demo5.py:24
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
void CrossCheckDiag(TMatrixDSym covMx, TMatrixD V, TMatrixD eigenvectors, TVectorD eigenvalues, double tolerance=1e-14)
std::vector< std::string > axis_labels
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:30
int nunivs
void FillFlavorContainers(std::string fileName, std::vector< std::string > comps, std::vector< std::string > comps2)
void plot_label(double x, double y, std::string label)
void SavePCAForCAFAna ( int  ncomponents,
double  inflate = 1. 
)

Definition at line 260 of file ppfx_syst_pca_fn.C.

References comps, genie::utils::style::Format(), hminus, hplus, MECModelEnuComparisons::i, nbins, pc_shifts_minus, pc_shifts_plus, and SplitHistograms().

Referenced by ppfx_syst_pca_fn().

261 {
262  TFile* saveHists = new TFile("ppfx_hadp_pc_shifts_fn_2018.root", "RECREATE");
263 
264  int ncomps = comps.size();
265  std::vector<std::string> histnames = {
266  "RHC_ND_numu",
267  "RHC_ND_nue",
268  "RHC_ND_anumu",
269  "RHC_ND_anue",
270  "FHC_ND_numu",
271  "FHC_ND_nue",
272  "FHC_ND_anumu",
273  "FHC_ND_anue",
274  "RHC_FD_numu",
275  "RHC_FD_nue",
276  "RHC_FD_anumu",
277  "RHC_FD_anue",
278  "FHC_FD_numu",
279  "FHC_FD_nue",
280  "FHC_FD_anumu",
281  "FHC_FD_anue",
282  };
283 
284  for(int pcIdx = 0; pcIdx < ncomponents; pcIdx++) {
285  std::vector<TH1D*> hplus = SplitHistograms(pc_shifts_plus[pcIdx], ncomps);
286  std::vector<TH1D*> hminus = SplitHistograms(pc_shifts_minus[pcIdx], ncomps);
287 
288  for(int hIdx = 0; hIdx < ncomps; hIdx++) {
289 
290  std::string pcIdxStr = TString::Format("%02d", pcIdx).Data();
291  std::string nplus = "ppfx_hadp_pc"+pcIdxStr+"_"+histnames[hIdx]+"_plussigma";
292  std::string nminus = "ppfx_hadp_pc"+pcIdxStr+"_"+histnames[hIdx]+"_minussigma";
293 
294  int nbins = hplus[hIdx]->GetNbinsX();
295  TH1D* ppfxplus = new TH1D(nplus.c_str(), nplus.c_str(), nbins, 0, 10);
296  TH1D* ppfxminus = new TH1D(nminus.c_str(), nminus.c_str(), nbins, 0, 10);
297 
298  for(int i = 1; i <= hplus[hIdx]->GetNbinsX(); i++){
299  ppfxplus->SetBinContent(i, inflate*(hplus[hIdx]->GetBinContent(i)-1.)+1.);
300  ppfxminus->SetBinContent(i, inflate*(hminus[hIdx]->GetBinContent(i)-1.)+1.);
301  }
302  ppfxplus->Write();
303  ppfxminus->Write();
304  }
305 
306  }
307  saveHists->cd();
308  saveHists->Close();
309 
310 }
std::vector< TH1D * > SplitHistograms(TH1D *hist, int nhists)
std::vector< TH1D * > pc_shifts_minus
TH1F * hplus
Definition: Xsec_final.C:75
std::vector< TH1D * > pc_shifts_plus
const int nbins
Definition: cellShifts.C:15
std::vector< std::string > comps
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
TH1F * hminus
Definition: Xsec_final.C:76
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
void SavePCAShifts ( int  ncomponents,
TH1D *  joint_hist_cv,
std::vector< TH1D * >  joint_hist_systs 
)

Definition at line 180 of file ppfx_syst_pca_fn.C.

References CombineHistograms(), comps, cv_areas, fn_cv, fn_shifts_plus, fn_shifts_ppfx, makeTrainCVSamples::int, pc_shifts_plus, Scale(), SplitHistograms(), art::to_string(), and ana::UniqueName().

Referenced by ppfx_syst_pca_fn().

182 {
183  TFile* saveHists = new TFile("check_cover_pc_fn_hadp_random.root", "RECREATE");
184 
185  TDirectory *cvdir = saveHists->mkdir("cv");
186  cvdir->cd();
187  TH1D* cv_save_hist = (TH1D*)joint_hist_cv->Clone("ppfx_cv");
188  cv_save_hist->Write();
189  TH1D* fn_cv_save = (TH1D*)fn_cv->Clone("fn_cv");
190  fn_cv_save->Write();
191  std::vector<TH1D*> hflavs = SplitHistograms(joint_hist_cv, comps.size());
192  for(int comp = 0; comp < (int)comps.size(); comp++)
193  hflavs[comp]->Scale(cv_areas[comp]);
194 
195  std::vector<TH1D*> near_flavs(hflavs.begin(), hflavs.begin()+8);
196  std::vector<TH1D*> far_flavs(hflavs.begin()+8, hflavs.end());
197 
198  TH1D* hnear = CombineHistograms(near_flavs);
199  hnear->SetName("near_cv");
200  TH1D* hfar = CombineHistograms(far_flavs);
201  hfar->SetName("far_cv");
202  hnear->Write();
203  hfar->Write();
204  saveHists->cd();
205 
206  int UnivIdx=0;
207  TDirectory *univdir = saveHists->mkdir("univs");
208  univdir->cd();
209  for(auto joint_hist_univ: joint_hist_systs) {
210  std::string name_univ = "ppfx_univ"+std::to_string(UnivIdx);
211  TH1D* joint_hist_syst_ratio = (TH1D*)joint_hist_univ->Clone(name_univ.c_str());
212  joint_hist_syst_ratio->Write();
213  UnivIdx++;
214  }
215  saveHists->cd();
216 
217  TDirectory *pcdir = saveHists->mkdir("shifts");
218  pcdir->cd();
219  int pcIdx = 0;
220  for(auto pc_shift:pc_shifts_plus) {
221  TH1D* ppfxplus = (TH1D*)pc_shift->Clone(UniqueName().c_str());
222  std::string nameplus_save = "ppfx"+std::to_string(pcIdx);
223  ppfxplus->SetName(nameplus_save.c_str());
224  ppfxplus->Write();
225 
226  if(pcIdx >= ncomponents) break;
227  pcIdx++;
228  }
229  saveHists->cd();
230 
231  TDirectory *fndir = saveHists->mkdir("fn_ratios");
232  fndir->cd();
233  int fnIdx = 0;
234  for(auto h_fn: fn_shifts_plus) {
235  TH1D* fn_ratio_save = (TH1D*) h_fn->Clone(UniqueName().c_str());
236  std::string name_fn = "fn_ratios_"+std::to_string(fnIdx);
237  fn_ratio_save->SetName(name_fn.c_str());
238  fn_ratio_save->Write();
239 
240  if(fnIdx >= ncomponents) break;
241  fnIdx++;
242  }
243  saveHists->cd();
244 
245  TDirectory *fnppfxdir = saveHists->mkdir("fn_ratios_ppfx");
246  fnppfxdir->cd();
247  int fnppfxIdx = 0;
248  for(auto h_fn: fn_shifts_ppfx) {
249  TH1D* fn_ratio_save = (TH1D*) h_fn->Clone(UniqueName().c_str());
250  std::string name_fn = "fn_ratios_ppfx_"+std::to_string(fnppfxIdx);
251  fn_ratio_save->SetName(name_fn.c_str());
252  fn_ratio_save->Write();
253  fnppfxIdx++;
254  }
255  saveHists->cd();
256  saveHists->Close();
257 
258 }
TH1D * fn_cv
std::vector< TH1D * > SplitHistograms(TH1D *hist, int nhists)
std::vector< TH1D * > fn_shifts_ppfx
TH1D * CombineHistograms(std::vector< TH1D * > hists)
std::vector< double > cv_areas
std::vector< TH1D * > pc_shifts_plus
std::vector< TH1D * > fn_shifts_plus
std::vector< std::string > comps
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
simulatedPeak Scale(1/simulationNormalisationFactor)
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:30

Variable Documentation

std::vector<std::string> comps
Initial value:
= {
"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",
}

Definition at line 27 of file ppfx_syst_pca_fn.C.

Referenced by ana::covmx::CovarianceMatrix::BuildFullCovMx(), skim::ParametersNumuCCpi::Compare(), skim::ParametersNue::Compare(), skim::ParametersNumu::Compare(), ana::CompareNDDataMCFromVector(), ConvertNormalBasis(), GetFNRatio(), plot_nd_spectra_2018(), PlotComp(), ppfx_make_systs(), ppfx_smooth_weights_make(), ppfx_smooth_weights_save(), ppfx_syst_pca_fn(), ana::covmx::CovarianceMatrix::Predict(), SavePCAForCAFAna(), SavePCAShifts(), and ana::covmx::CovarianceMatrix::SetNUniverses().

std::vector<std::string> comps2
Initial value:
= {
"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",
}

Definition at line 45 of file ppfx_syst_pca_fn.C.

Referenced by ppfx_syst_pca_fn().

std::vector<double> cv_areas
std::vector<double> cv_areas_fn

Definition at line 74 of file ppfx_syst_pca_fn.C.

Referenced by ConvertFNBasis(), and ConvertNormalBasis().

std::vector<TH1D*> cv_comps

Definition at line 69 of file ppfx_syst_pca_fn.C.

Referenced by ConvertFNBasis(), and FillFlavorContainers().

std::vector<TH1D*> cv_comps_fn

Definition at line 73 of file ppfx_syst_pca_fn.C.

Referenced by ConvertFNBasis(), and ppfx_syst_pca_fn().

TH1D* fn_cv

Definition at line 83 of file ppfx_syst_pca_fn.C.

Referenced by ppfx_check_var(), ppfx_syst_pca_fn(), and SavePCAShifts().

std::vector<TH1D*> fn_shifts_minus

Definition at line 86 of file ppfx_syst_pca_fn.C.

Referenced by ppfx_syst_pca_fn().

std::vector<TH1D*> fn_shifts_plus

Definition at line 85 of file ppfx_syst_pca_fn.C.

Referenced by ppfx_syst_pca_fn(), and SavePCAShifts().

std::vector<TH1D*> fn_shifts_ppfx

Definition at line 84 of file ppfx_syst_pca_fn.C.

Referenced by ppfx_syst_pca_fn(), and SavePCAShifts().

int nunivs = 100

Definition at line 66 of file ppfx_syst_pca_fn.C.

Referenced by ConvertFNBasis(), FillFlavorContainers(), and ppfx_syst_pca_fn().

std::vector<TH1D*> pc_shifts_minus

Definition at line 79 of file ppfx_syst_pca_fn.C.

Referenced by FillPCAContainers(), ppfx_syst_pca_fn(), and SavePCAForCAFAna().

std::vector<TH1D*> pc_shifts_plus
int rebin_f = 10

Definition at line 65 of file ppfx_syst_pca_fn.C.

std::vector<std::vector<TH1D*> > univs_comps

Definition at line 71 of file ppfx_syst_pca_fn.C.

Referenced by ConvertFNBasis(), and FillFlavorContainers().

std::vector<std::vector<TH1D*> > univs_comps_fn

Definition at line 75 of file ppfx_syst_pca_fn.C.

Referenced by ConvertFNBasis(), and ppfx_syst_pca_fn().