MatrixValidation.C
Go to the documentation of this file.
1 /// MatrixValidation.C by J. Hewes <jhewes15@fnal.gov>
2 
3 #ifdef __CINT__
4 
5 #include <iostream>
6 
7 void MatrixValidation(bool rhc = false) {
8  throw std::runtime_error("You must run in compiled mode.");
9 }
10 
11 #else
12 
13 #include "TLine.h"
14 
19 
21 
22 using namespace ana;
23 
24 TLatex* MiscText(double x, double y, double size, TString text);
25 std::string GetName(int i);
26 
27 void MatrixValidation(bool rhc = false) {
28 
29  // ROOT guard
30  DontAddDirectory guard;
31 
32  // Preamble
33  FormatCanvas();
34 
36  double fd_pot = rhc ? kAna2018RHCPOT : kAna2018FHCPOT;
37  double fd_livetime = rhc ? kAna2018RHCLivetime : kAna2018FHCLivetime;
38 
39  std::string polarity = rhc ? "RHC" : "FHC";
40 
42  SetAngles(calc);
43  osc::OscCalcSterileTrivial* calc_trivial
44  = new osc::OscCalcSterileTrivial(*calc);
45 
46  // Load and concatenate predictions
47  PredictionConcat* sim = GetDefaultSimulation(rhc);
48 
49  // Load all the covariance matrices
50  std::map<const ISyst*, CovarianceMatrix*> m;
51  for (const ISyst* syst : GetNus18Systs(false, "base", "xsec")) {
52  CovarianceMatrix* gen = LoadCovMx("covmx_fhc_xsec", "CovMx_Both_" + syst->ShortName());
53  m.insert(std::pair<const ISyst*, CovarianceMatrix*>(syst, gen));
54  }
55  for (const ISyst* syst : GetNus18Systs(false, "base", "nonxsec")) {
56  CovarianceMatrix* gen = LoadCovMx("covmx_fhc_nonxsec", "CovMx_Both_" + syst->ShortName());
57  m.insert(std::pair<const ISyst*, CovarianceMatrix*>(syst, gen));
58  }
59  CovarianceMatrix* gen = LoadCovMx("covmx_fhc_ppfx", "CovMx_Both_PPFX");
60  m.insert(std::pair<const ISyst*, CovarianceMatrix*>(nullptr, gen));
61  CovarianceMatrix* matrix_combined = new CovarianceMatrix(sim, calc_trivial, {}, fd_pot, 10);
62  CombineMatrices(matrix_combined, m, sim, calc);
63 
64 
65  // Now draw these matrices
66  TCanvas* c = new TCanvas();
67  std::string filename = "plots/covmx/" + polarity + "_Plots.root";
68  TFile* plotfile = new TFile(filename.c_str(), "recreate");
69 
70  for (auto matrix : m) {
71 
72  // Full matrix
73  TH2D* h_full = matrix.second->GetFullCovMxTH2();
74  int nbins = h_full->GetNbinsX();
75  double last_edge = sim->GetBinEdges().back();
76  double max = 13. * last_edge;
77  h_full->Draw("colz");
78 
79  for (int i = 1; i < 13; ++i) {
80  double pos = i * last_edge;
81  TLine* lx = new TLine(0, pos, max, pos);
82  lx->Draw();
83  TLine* ly = new TLine(pos, 0, pos, max);
84  ly->Draw();
85  }
86 
87  for (int i = 0; i < 13; ++i) {
88  double pos = 0.16 + (0.72*(i/13.));
89  auto l1 = MiscText(pos,0.88,0.04,GetName(i).c_str());
90  l1->SetTextAngle(38);
91  auto l2 = MiscText(0.12,pos,0.04,GetName(i).c_str());
92  l2->SetTextAlign(31);
93  }
94 
95  h_full->SetNdivisions(0,"X");
96  h_full->SetNdivisions(0,"Y");
97 
98  std::string syst_name = (matrix.first == nullptr) ? "PPFX" : matrix.first->ShortName();
99  std::string plot_name = "CovMx_" + polarity + "_" + syst_name + "_Full";
100  std::string plot_path = "plots/covmx/" + plot_name + ".png";
101  c->Write(plot_name.c_str());
102  c->SaveAs(plot_path.c_str());
103  c->Clear();
104 
105  // Oscillated matrix
106  TH1D* nominal = sim->Predict(calc).ToTH1(fd_pot);
107  nbins = nominal->GetNbinsX();
108  TVectorD v_pred(nbins);
109  for (int i = 0; i < nbins; ++i)
110  v_pred(i) = nominal->GetBinContent(i+1);
111  TH2D* h_osc = matrix.second->GetCovMxRelativeTH2(v_pred);
112  h_osc->Draw("colz");
113  h_osc->GetYaxis()->SetTitle("");
114  h_osc->GetXaxis()->SetTitle("Deposited energy bins");
115  h_osc->GetXaxis()->CenterTitle();
116 
117  Preliminary();
118 
119  double max_e = 20;
120 
121  TLine lx(0,max_e,2*max_e,max_e);
122  lx.Draw();
123  TLine ly(max_e,0,max_e,2*max_e);
124  ly.Draw();
125  MiscText(0.32,0.88,0.06,"ND");
126  MiscText(0.68,0.88,0.06,"FD");
127  MiscText(0.04,0.3,0.06,"ND");
128  MiscText(0.04,0.7,0.06,"FD");
129  plot_name = "CovMx_" + polarity + "_" + syst_name + "_Oscillated";
130  plot_path = "plots/covmx/" + plot_name + ".png";
131  c->Write(plot_name.c_str());
132  c->SaveAs(plot_path.c_str());
133  c->Clear();
134  }
135 
136  //
137  // EFFECTIVE PULL EXTRACTION
138  //
139 
140  // Get vector of all ISysts
141  std::vector<const ISyst*> systs;
142  std::vector<std::string> syst_types = { "xsec", "nonxsec" };
143  for (std::string t : syst_types) {
144  std::vector<const ISyst*> s = GetNus18Systs(rhc, "base", t);
145  systs.insert(systs.end(), s.begin(), s.end());
146  }
147 
148  // Loop over this vector
149  for (const ISyst* syst : systs) {
150 
151  // Generate new empty matrix
152  CovarianceMatrix* mx = new CovarianceMatrix(sim, calc_trivial, {}, fd_pot, 10);
153 
154  for (auto matrix : m) {
155 
156  std::string syst_name = matrix.first == nullptr ? "PPFX" : matrix.first->ShortName();
157 
158  // Drop this systematic out
159  if (syst->ShortName() == syst_name) continue;
160 
161  // Add the other systematics to the new matrix
162  mx->AddMatrix(matrix.second);
163  }
164 
165  // Predict matrix to ensure we've reset it correctly
166  mx->PredictCovMx(sim, calc);
167 
168  // Draw matrix
169  auto nominal = sim->Predict(calc).ToTH1(fd_pot);
170  double nbins = nominal->GetNbinsX();
171  TVectorD v_pred(nbins);
172  for (int i = 0; i < nbins; ++i)
173  v_pred(i) = nominal->GetBinContent(i+1);
174  TH2D* h_osc_rel = mx->GetCovMxRelativeTH2(v_pred);
175  h_osc_rel->Draw("colz");
176  h_osc_rel->GetYaxis()->SetTitle("");
177  h_osc_rel->GetXaxis()->SetTitle("Deposited energy bins");
178  h_osc_rel->GetXaxis()->CenterTitle();
179 
180  // Preliminary();
181 
182  double max_e = 20;
183 
184  TLine lx(0, max_e, 2 * max_e, max_e);
185  lx.Draw();
186  TLine ly(max_e, 0, max_e, 2 * max_e);
187  ly.Draw();
188  MiscText(0.32,0.88,0.06,"ND");
189  MiscText(0.68,0.88,0.06,"FD");
190  MiscText(0.04,0.3,0.06,"ND");
191  MiscText(0.04,0.7,0.06,"FD");
192  c->SaveAs(Form("plots/covmx/CovMxWithout%sOscRel.png",syst->ShortName().c_str()));
193  c->SaveAs(Form("plots/covmx/CovMxWithout%sOscRel.pdf",syst->ShortName().c_str()));
194  c->Clear();
195  delete h_osc_rel;
196 
197  // Now extract the pull manually
198  SystShifts shift;
199  double pull = -5;
200  for (int i = 0; i < 21; ++i) {
201  shift.SetShift(syst, pull);
202  pull += 0.5;
203  }
204  }
205 }
206 
207 //------------------------------------------------------------------------------
208 TLatex* MiscText(double x, double y, double size, TString text)
209 {
210  TLatex *l = new TLatex(x, y, text);
211  l->SetNDC();
212  l->SetTextSize(size);
213  l->SetTextColor(kBlack);
214  l->Draw();
215  return l;
216 }
217 
218 //------------------------------------------------------------------------------
220  if (i==0) return std::string("#nu_{e} #rightarrow #nu_{e}");
221  else if (i==1) return std::string("#bar{#nu}_{e} #rightarrow #bar{#nu}_{e}");
222  else if (i==2) return std::string("#nu_{e} #rightarrow #nu_{#mu}");
223  else if (i==3) return std::string("#bar{#nu}_{e} #rightarrow #bar{#nu}_{#mu}");
224  else if (i==4) return std::string("#nu_{e} #rightarrow #nu_{#tau}");
225  else if (i==5) return std::string("#bar{#nu}_{e} #rightarrow #bar{#nu}_{#tau}");
226  else if (i==6) return std::string("#nu_{#mu} #rightarrow #nu_{e}");
227  else if (i==7) return std::string("#bar{#nu}_{#mu} #rightarrow #bar{#nu}_{e}");
228  else if (i==8) return std::string("#nu_{#mu} #rightarrow #nu_{#mu}");
229  else if (i==9) return std::string("#bar{#nu}_{#mu} #rightarrow #bar{#nu}_{#mu}");
230  else if (i==10) return std::string("#nu_{#mu} #rightarrow #nu_{#tau}");
231  else if (i==11) return std::string("#bar{#nu}_{#mu} #rightarrow #bar{#nu}_{#tau}");
232  else if (i==12) return std::string("NC");
233  std::cerr << "What are you doing???" << std::endl;
234  return "no";
235 }
236 
237 #endif
238 
void MatrixValidation(bool rhc=false)
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
Adapt the PMNS_Sterile calculator to standard interface.
OStream cerr
Definition: OStream.cxx:7
Version of OscCalcSterile that always returns probability of 1.
string filename
Definition: shutoffs.py:106
const double kAna2018SensitivityFHCNDPOT
Definition: Exposures.h:210
osc::OscCalcDumb calc
void FormatCanvas()
Canvas formatting utility.
void SetAngles(osc::OscCalcSterile *calc)
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
const double kAna2018RHCPOT
Definition: Exposures.h:208
const XML_Char * s
Definition: expat.h:262
CovarianceMatrix * LoadCovMx(std::string file, std::string name)
Function to load covariance matrix.
const int nbins
Definition: cellShifts.C:15
osc::OscCalcSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
Definition: Calcs.cxx:97
std::string GetName(int i)
TLatex * MiscText(double x, double y, double size, TString text)
Example routine for calculating the "ultimate e-m mother" of a particle in a simulated event...
Definition: FillTruth.h:16
void Preliminary()
const double kAna2018SensitivityRHCNDPOT
Definition: Exposures.h:211
const double kAna2018FHCPOT
Definition: Exposures.h:207
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:80
std::vector< const ISyst * > GetNus18Systs(bool rhc, std::string det_type, std::string syst_type)
Definition: Nus18Systs.cxx:439
const double kAna2018FHCLivetime
Definition: Exposures.h:213
const double kAna2018RHCLivetime
Definition: Exposures.h:214
void CombineMatrices(CovarianceMatrix *gen, std::map< const ISyst *, CovarianceMatrix * > m, std::vector< IPrediction * > preds, std::vector< double > pot, osc::IOscCalcAdjustable *calc)
Function to combine covariance matrices.
enum BeamMode string