13 #include "TMatrixDSymEigen.h" 14 #include "TMatrixDSym.h" 27 std::vector<std::string>
comps = {
30 "nd_RHC_trueE_anti_numu",
31 "nd_RHC_trueE_anti_nue",
34 "nd_FHC_trueE_anti_numu",
35 "nd_FHC_trueE_anti_nue",
36 "fd_RHC_trueE_numu_ns",
37 "fd_RHC_trueE_nue_ns",
38 "fd_RHC_trueE_anti_numu_ns",
39 "fd_RHC_trueE_anti_nue_ns",
40 "fd_FHC_trueE_numu_ns",
41 "fd_FHC_trueE_nue_ns",
42 "fd_FHC_trueE_anti_numu_ns",
43 "fd_FHC_trueE_anti_nue_ns",
46 "RHC_univ_nd_trueE_numu",
47 "RHC_univ_nd_trueE_nue",
48 "RHC_univ_nd_trueE_anti_numu",
49 "RHC_univ_nd_trueE_anti_nue",
50 "FHC_univ_nd_trueE_numu",
51 "FHC_univ_nd_trueE_nue",
52 "FHC_univ_nd_trueE_anti_numu",
53 "FHC_univ_nd_trueE_anti_nue",
54 "RHC_univ_fd_trueE_numu_ns",
55 "RHC_univ_fd_trueE_nue_ns",
56 "RHC_univ_fd_trueE_anti_numu_ns",
57 "RHC_univ_fd_trueE_anti_nue_ns",
58 "FHC_univ_fd_trueE_numu_ns",
59 "FHC_univ_fd_trueE_nue_ns",
60 "FHC_univ_fd_trueE_anti_numu_ns",
61 "FHC_univ_fd_trueE_anti_nue_ns",
90 TFile* fload =
new TFile(fileName.c_str(),
"READ");
92 for(
int compIdx = 0; compIdx < (
int)comps.size(); ++compIdx){
93 const char* comp = comps[compIdx].c_str();
97 TH1D* temp_hist = (TH1D*)spec->ToTH1(spec->POT());
99 cv_areas.push_back(temp_hist->Integral());
100 temp_hist->Scale(1./temp_hist->Integral());
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();
111 TH1D* hist_syst = (TH1D*)spec_syst->ToTH1(spec_syst->POT());
113 hist_syst->Scale(1./
cv_areas[compIdx]);
114 temp_hist_systs.push_back(hist_syst);
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());
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());
131 normal_basis.push_back(h_normal);
133 for(
int fnIdx = 0; fnIdx < (
int)fn_hists.size(); fnIdx++){
136 far_h->Multiply(normal_basis[fnIdx]);
138 normal_basis.push_back(far_h);
145 TH1D* joint_hist_cv_fn, std::vector<TH1D*> joint_hist_systs_fn,
146 TH1D* joint_hist_cv, std::vector<TH1D*> joint_hist_systs)
154 TH1D* retplus_fn =
new TH1D(x);
156 retplus_fn->Add(joint_hist_cv_fn);
159 retplus->SetName(nameplus.c_str());
160 retplus->Divide(joint_hist_cv);
161 retplus->GetYaxis()->SetRangeUser(0.7, 1.3);
166 TH1D* retminus_fn =
new TH1D(x);
168 retminus_fn->Add(joint_hist_cv_fn);
171 retminus->SetName(nameminus.c_str());
172 retminus->Divide(joint_hist_cv);
173 retminus->GetYaxis()->SetRangeUser(0.7, 1.3);
181 TH1D* joint_hist_cv, std::vector<TH1D*> joint_hist_systs)
183 TFile* saveHists =
new TFile(
"check_cover_pc_fn_hadp_random.root",
"RECREATE");
185 TDirectory *cvdir = saveHists->mkdir(
"cv");
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");
192 for(
int comp = 0; comp < (
int)
comps.size(); comp++)
195 std::vector<TH1D*> near_flavs(hflavs.begin(), hflavs.begin()+8);
196 std::vector<TH1D*> far_flavs(hflavs.begin()+8, hflavs.end());
199 hnear->SetName(
"near_cv");
201 hfar->SetName(
"far_cv");
207 TDirectory *univdir = saveHists->mkdir(
"univs");
209 for(
auto joint_hist_univ: joint_hist_systs) {
211 TH1D* joint_hist_syst_ratio = (TH1D*)joint_hist_univ->Clone(name_univ.c_str());
212 joint_hist_syst_ratio->Write();
217 TDirectory *pcdir = saveHists->mkdir(
"shifts");
221 TH1D* ppfxplus = (TH1D*)pc_shift->Clone(
UniqueName().c_str());
223 ppfxplus->SetName(nameplus_save.c_str());
226 if(pcIdx >= ncomponents)
break;
231 TDirectory *fndir = saveHists->mkdir(
"fn_ratios");
235 TH1D* fn_ratio_save = (TH1D*) h_fn->Clone(
UniqueName().c_str());
237 fn_ratio_save->SetName(name_fn.c_str());
238 fn_ratio_save->Write();
240 if(fnIdx >= ncomponents)
break;
245 TDirectory *fnppfxdir = saveHists->mkdir(
"fn_ratios_ppfx");
249 TH1D* fn_ratio_save = (TH1D*) h_fn->Clone(
UniqueName().c_str());
251 fn_ratio_save->SetName(name_fn.c_str());
252 fn_ratio_save->Write();
262 TFile* saveHists =
new TFile(
"ppfx_hadp_pc_shifts_fn_2018.root",
"RECREATE");
264 int ncomps =
comps.size();
265 std::vector<std::string> histnames = {
284 for(
int pcIdx = 0; pcIdx < ncomponents; pcIdx++) {
288 for(
int hIdx = 0; hIdx < ncomps; hIdx++) {
291 std::string nplus =
"ppfx_hadp_pc"+pcIdxStr+
"_"+histnames[hIdx]+
"_plussigma";
292 std::string nminus =
"ppfx_hadp_pc"+pcIdxStr+
"_"+histnames[hIdx]+
"_minussigma";
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);
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.);
316 for(
int comp = 0; comp < ncomps; comp++)
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++){
324 hfn->Divide(near_flavs[flavIdx]);
325 fn_flavs.push_back(hfn);
343 TH1D*
FNBasisHelper(TH1D* near, TH1D* far,
double near_norm,
double far_norm)
345 TH1D* fn = (TH1D*)far->Clone(
UniqueName().c_str());
347 fn->Scale(far_norm/near_norm);
354 for(
int nearIdx = 0; nearIdx < 8; nearIdx++) {
356 near->Scale(nonfn_scale);
360 for(
int farIdx = 8; farIdx < 16; farIdx++) {
363 fn->Scale(1./fn->Integral());
367 for(
int UnivIdx = 0; UnivIdx < (
int)
nunivs; UnivIdx++){
368 std::vector<TH1D*> temp_fn_univ;
369 for(
int nearIdx = 0; nearIdx < 8; nearIdx++) {
371 near_univ->Scale(nonfn_scale);
372 temp_fn_univ.push_back(near_univ);
374 for(
int farIdx = 8; farIdx < 16; farIdx++) {
377 temp_fn_univ.push_back(fn_univ);
396 std::vector <TH1D*> joint_hist_systs_fn;
399 joint_hist_systs_fn.push_back(joint_hist_univ);
404 int joint_nbins = joint_hist_cv_fn->GetNbinsX();
407 const std::vector<TH1D*> ppfx(joint_hist_systs_fn.begin(), joint_hist_systs_fn.end());
410 std::unique_ptr<TMatrixDSym> covMx =
GetCovMx(ppfx, joint_hist_cv_fn);
413 TMatrixDSymEigen cov(*covMx.get());
416 eigenvectors.Transpose(V);
417 TVectorD eigenvalues = cov.GetEigenValues();
423 std::vector<TH1D*> joint_hist_systs;
424 std::for_each(joint_hist_systs_fn.begin(), joint_hist_systs_fn.end(), [&](TH1D* joint_h){
429 FillPCAContainers(eigenvectors, eigenvalues, joint_hist_cv_fn, joint_hist_systs_fn, joint_hist_cv, joint_hist_systs);
432 int ncomponents = 100;
439 for(
int i = 0; i < ncomponents; i++){
441 fn_plus->Multiply(joint_hist_cv);
445 fn_minus->Multiply(joint_hist_cv);
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",
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",
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",
473 std::vector <plot_h> h_cv_systs;
474 std::vector <std::pair<plot_h, TH1*>> h_cv_systs_ratios;
476 for(
auto hist: joint_hist_systs){
478 h_cv_systs.push_back(univ);
482 plot_h cv = {joint_hist_cv,
"hist same", kBlack,
false,
"",
"l"};
483 h_cv_systs.push_back(cv);
486 TCanvas *
c1 =
new TCanvas(
"c1",
"c1", 200, 10, 800, 600);
493 c1->Print(
"ppfx_diag_plots/ppfx_hadp_spectra.pdf");
495 TCanvas *
c2 =
new TCanvas(
"c2",
"c2", 200, 10, 800, 600);
500 draw_vlines(
comps.size(), joint_hist_cv->GetXaxis()->GetBinLowEdge(joint_nbins+1), 0, 2);
502 c2->Print(
"ppfx_diag_plots/ppfx_hadp_spectra_ratios.pdf");
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]);
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]);
518 gStyle->SetPalette(55);
519 TCanvas *
c3 =
new TCanvas(
"c3",
"c3", 200, 10, 800, 600);
522 histCov->Draw(
"COLZ");
524 c3->SetRightMargin(0.125);
527 c3->Print(
"ppfx_diag_plots/ppfx_hadp_covariance.pdf");
529 TCanvas *c3_fn =
new TCanvas(
"c3_fn",
"c3_fn", 200, 10, 800, 600);
532 histCov_fn->Draw(
"COLZ");
534 c3_fn->SetRightMargin(0.125);
537 c3_fn->Print(
"ppfx_diag_plots/ppfx_hadp_covariance_fn.pdf");
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);
546 histEigen->Draw(
"l");
549 c4->Print(
"ppfx_diag_plots/ppfx_hadp_eigenvalues.pdf");
551 int nppfx_diag_plots = 5;
555 std::vector<plot_h> far_near_shifts = {
556 {
fn_cv,
"hist same", kBlack,
true,
"CV",
"l"},
561 std::vector<std::pair<plot_h, TH1*>> double_ratios;
562 std::for_each(far_near_shifts.begin(), far_near_shifts.end(),
566 std::vector<plot_h> pc_shifts;
574 TCanvas *c5 =
new TCanvas(
"c5",
"c5", 200, 10, 800, 600);
578 plot_label(0.5, 0.93, Form(
"Proportion of Variance : %lg", eigenvalues[
component]/eigenvalues.Sum()));
580 draw_vlines(
comps.size(), joint_hist_cv->GetXaxis()->GetBinLowEdge(joint_nbins+1), 0.7, 1.3);
583 c5->Print(filename1.c_str());
585 TCanvas *c6 =
new TCanvas(
"c6",
"c6", 200, 10, 800, 600);
587 TPad *padUp =
new TPad(
"padUp",
"padUp",0.0,0.30,1.,1.00);
590 padUp->SetBottomMargin(0.00);
591 padUp->SetBorderSize(0);
595 plot_hists(far_near_shifts, fn_title,
false,
false);
604 TPad *padDn =
new TPad(
"padDn",
"padDn",0.0,0.0,1.,0.30);
607 padDn->SetBottomMargin(0.30);
608 padDn->SetTopMargin(0.00);
609 padDn->SetBorderSize(0);
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.);
623 c6->Print(filename2.c_str());
TH1D * FNBasisHelper(TH1D *near, TH1D *far, double near_norm, double far_norm)
std::vector< std::string > comps2
TH1D * GetTotFNRatio(TH1D *joint_hist)
std::vector< TH1D * > SplitHistograms(TH1D *hist, int nhists)
std::vector< std::vector< int > > rebin_f
Cuts and Vars for the 2020 FD DiF Study.
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< double > cv_areas
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)
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
std::vector< TH1D * > cv_comps
std::vector< TH1D * > fn_shifts_plus
TH1D * ConvertNormalBasis(TH1D *joint_hist_fn)
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
std::vector< std::vector< TH1D * > > univs_comps
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
void ConvertFNBasis(double nonfn_scale=1e-6)
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)
std::vector< double > cv_areas_fn
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)
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
std::string to_string(ModuleType mt)
void CrossCheckDiag(TMatrixDSym covMx, TMatrixD V, TMatrixD eigenvectors, TVectorD eigenvalues, double tolerance=1e-14)
std::vector< std::string > axis_labels
simulatedPeak Scale(1/simulationNormalisationFactor)
std::string UniqueName()
Return a different string each time, for creating histograms.
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)