PlotCovariance.C
Go to the documentation of this file.
1 /// \title PlotCovariance.C
2 /// \author Connor Johnson
3 /// \date Jan 16, 2020
4 /// \brief A script which takes covariance matrices, and constructs our final pretty plots for NumuCC
5 
6 #include <iostream>
7 #include <string>
8 #include <algorithm>
9 
10 // ROOT includes
11 #include "TFile.h"
12 #include "TH1.h"
13 #include "TH2.h"
14 #include "TRatioPlot.h"
15 #include "THStack.h"
16 #include "TCanvas.h"
17 #include "TString.h"
18 #include "TSpline.h"
19 #include "TLegend.h"
20 #include "TLegendEntry.h"
21 
22 #include "TMatrixD.h"
23 #include "TArray.h"
24 #include "TStyle.h"
25 
26 #include "TMinuit.h"
27 
28 // Needed ONLY for legend placement
29 #include "CAFAna/Analysis/Plots.h"
30 
31 string cv_name = "cv";
32 
33 // Binning logic
34 
35 map<int,int> MuTperCosBins =
36 {
37  {1, 6 },
38  {2, 7 },
39  {3, 8 },
40  {4, 9 },
41  {5, 9 },
42  {6, 9 },
43  {7, 13},
44  {8, 14},
45  {9, 17},
46  {10, 20},
47  {11, 20},
48  {12, 20},
49  {13, 20}
50 };
51 
52 const std::vector<double> angBinEdges {
53  -1.0, 0.50, 0.56, 0.62, 0.68,
54  0.74, 0.80, 0.85, 0.88, 0.91,
55  0.94, 0.96, 0.98, 0.99, 1.0
56 };
57 
58 const std::vector<double> angBinEdges_reduced {
59  0.50, 0.56, 0.62, 0.68,
60  0.74, 0.80, 0.85, 0.88, 0.91,
61  0.94, 0.96, 0.98, 0.99, 1.0
62 };
63 
64 const std::vector<double> mukeBinEdges {
65  -10.0,
66  0.5, 0.6, 0.7, 0.8, 0.9,
67  1.0, 1.1, 1.2, 1.3, 1.4,
68  1.5, 1.6, 1.7, 1.8, 1.9,
69  2.0, 2.1, 2.2, 2.3, 2.4,
70  2.5, 120
71 };
72 
73 const std::vector<double> mukeBinEdges_reduced {
74  0.5, 0.6, 0.7, 0.8, 0.9,
75  1.0, 1.1, 1.2, 1.3, 1.4,
76  1.5, 1.6, 1.7, 1.8, 1.9,
77  2.0, 2.1, 2.2, 2.3, 2.4,
78  2.5
79 };
80 
81 const std::vector<double> enuBinEdges{
82  // -10,
83  //0.0, 0.25,
84  0.50, 0.75,
85  1.0, 1.25, 1.50, 1.75,
86  2.0, 2.25, 2.50, 2.75,
87  3.0, 3.25, 3.50, 3.75,
88  4.0, 4.25, 4.50, 4.75,
89  5.0//, 120
90 };
91 
92 const std::vector<double> q2BinEdges {
93  0.00,
94  0.05, 0.10, 0.15, 0.20, 0.25,
95  0.30, 0.40, 0.50, 0.60, 0.75,
96  0.90, 1.10, 1.40, 1.80, 2.8//,
97  //120.
98 };
99 
100 // Axis titles
101 string mukeTitle = "T_{#mu} (GeV)";
102 string muAngTitle = "cos(#theta_{#mu})";
103 string enuTitle = "E_{#nu} (GeV)";
104 string q2Title = "q^{2} (GeV^{2})";
105 
106 // XSec titles
107 string xsec2DTitle = "#frac{d^{2} #sigma}{d cos(#theta_{#mu}) d T_{#mu}} (cm^{2} / GeV / nucleon #times 10^{39})";
108 string xsecENuTitle = "#sigma / E_{#nu} (cm^{2} / GeV / nucleon #times 10^{39})";
109 string xsecQ2Title = "#frac{d #sigma}{d q^{2}} (cm^{2} / GeV^{2} / nucleon #times 10^{39})";
110 
111 TH2* Construct2D(TH1* histIn)
112 {
113  TH2* histOut = new TH2D((histIn->GetName() + string("_2D")).c_str(), "", angBinEdges_reduced.size() - 1, angBinEdges_reduced.data(), mukeBinEdges_reduced.size() - 1, mukeBinEdges_reduced.data());
114 
115  int counter = 1;
116  int NbinX = histOut->GetXaxis()->GetNbins();
117  int NbinY = histOut->GetYaxis()->GetNbins();
118 
119  // Logic reverse-engineered from covariance making script
120  for(int ibinx = 1; ibinx <= NbinX; ++ibinx){ // Throw out bin [-1, 0.5] in costheta
121  for(int ibiny = 1; ibiny <= NbinY; ++ibiny){ // Throw out first _and_ last bin of MuKE (-10, 0.5), (2.5, 120)
122  // Skip once we've reached how many MuKE bins to include in each costheta bin
123  if(MuTperCosBins.find(ibinx) == MuTperCosBins.end() || ibiny > MuTperCosBins[ibinx])
124  continue;
125  histOut->SetBinContent(ibinx, ibiny, histIn->GetBinContent(counter));
126  histOut->SetBinError(ibinx, ibiny, histIn->GetBinError(counter));
127  counter++;
128  }
129  }
130  histOut->GetXaxis()->SetRangeUser(0.5, 1.0);
131  histOut->GetYaxis()->SetRangeUser(0.5, 2.5);
132  return histOut;
133 }
134 
135 TH1* ConstructENu(TH1* histIn)
136 {
137  assert(histIn->GetNbinsX() == (int) enuBinEdges.size() - 1);
138  TH1 * histOut = new TH1D((histIn->GetName() + string("_enuaxis")).c_str(), "", enuBinEdges.size() - 1, enuBinEdges.data());
139  for (int ibinx = 1; ibinx <= histOut->GetNbinsX(); ++ibinx){
140  histOut->SetBinContent(ibinx, histIn->GetBinContent(ibinx));
141  histOut->SetBinError(ibinx, histIn->GetBinError(ibinx));
142  }
143  return histOut;
144 }
145 
146 TH1* ConstructQ2(TH1* histIn)
147 {
148  assert(histIn->GetNbinsX() == (int) q2BinEdges.size() - 1);
149  TH1 * histOut = new TH1D((histIn->GetName() + string("_q2axis")).c_str(), "", q2BinEdges.size() - 1, q2BinEdges.data());
150  for (int ibinx = 1; ibinx <= histOut->GetNbinsX(); ++ibinx){
151  histOut->SetBinContent(ibinx, histIn->GetBinContent(ibinx));
152  histOut->SetBinError(ibinx, histIn->GetBinError(ibinx));
153  }
154  return histOut;
155 }
156 
157 // Function to create the map from covariance bin # to global bin number in 2D
158 TH1I * MakeCovTo2DBin(TH2 * hist, int nbins)
159 {
160  TH1I * arr = new TH1I("covTo2DBin", "", nbins, 1, nbins + 1);
161  int bincounter = 1;
162  int NbinX = angBinEdges_reduced.size();
163  int NbinY = mukeBinEdges_reduced.size();
164  for(int ibinx = 1; ibinx <= NbinX; ++ibinx){ // Should this not be 2?
165  for(int ibiny = 1; ibiny <= NbinY; ++ibiny){
166  if(MuTperCosBins.find(ibinx) == MuTperCosBins.end() || ibiny > MuTperCosBins[ibinx])
167  continue;
168  arr->SetBinContent(bincounter, hist->GetBin(ibinx, ibiny));
169  bincounter++;
170  }
171  }
172  --bincounter; // Fix an off-by-one error for the very last valid bin incrementing.
173  cout << "Counter=" << bincounter << " NBins=" << nbins << endl;
174  assert(bincounter == nbins);
175  return arr;
176 }
177 
178 // Function to map x and y bin number to covariance bin number.
180 {
181  TH2I * binmap = new TH2I("hist2DBinToCov", "",
182  angBinEdges_reduced.size() - 1, angBinEdges_reduced.data(),
183  mukeBinEdges_reduced.size() - 1, mukeBinEdges_reduced.data());
184  int bincounter = 1;
185  int NbinX = angBinEdges_reduced.size();
186  int NbinY = mukeBinEdges_reduced.size();
187  for(int ibinx = 1; ibinx <= NbinX; ++ibinx){ // Should this not be 2?
188  for(int ibiny = 1; ibiny <= NbinY; ++ibiny){
189  if(MuTperCosBins.find(ibinx) == MuTperCosBins.end() || ibiny > MuTperCosBins[ibinx])
190  continue;
191  binmap->SetBinContent(ibinx, ibiny, bincounter);
192  bincounter++;
193  }
194  }
195  --bincounter; // Fix an off-by-one error for the very last valid bin incrementing.
196  cout << "Counter=" << bincounter << endl;
197  return binmap;
198 }
199 
200 // Function to go from global bin number to bin in the covariance matrix
201 int GetCovBin(TH1I * covTo2DBin, int globalBin)
202 {
203  assert(covTo2DBin && covTo2DBin->GetNbinsX() > 0);
204  for(int i = 1; i <= covTo2DBin->GetNbinsX(); ++i)
205  if (covTo2DBin->GetBinContent(i) == globalBin)
206  return i;
207  return -1;
208 }
209 
210 // Set the errors of a histogram to the sqrt of diagonals of the corresponding covariance matrix
211 TH1 * SetErrorFromCovDiag(TH1* hist, TH2* mat)
212 {
213  assert(hist->GetNbinsX() == mat->GetNbinsX());
214  for (int ibin = 1; ibin <= mat->GetNbinsX(); ++ibin)
215  {
216  hist->SetBinError(ibin, sqrt(mat->GetBinContent(ibin, ibin))); // Fill with sqrt diagonals
217  }
218  return hist;
219 }
220 
221 // Set the errors of a histogram to the sqrt of diagonals of the corresponding covariance matrix
223 {
224  assert(hist->GetNbinsX() == mat.GetNrows());
225  for (int irow = 0; irow < mat.GetNrows(); ++irow)
226  {
227  hist->SetBinError(irow+1, sqrt(mat[irow][irow])); // Fill with sqrt diagonals
228  }
229  return hist;
230 }
231 
232 // Helpers for getting
233 enum EBinEdge{
236 };
237 
238 int GetBinByEdge(TAxis* axis, double binEdge, EBinEdge findBinBy = kBinLowEdge){
239  double epsilon = 1E-5;
240  int foundBin = -1;
241  for(int ibin = 1; ibin <= axis->GetNbins(); ++ibin){
242  double ibinEdge;
243  if (findBinBy == kBinLowEdge)
244  ibinEdge = axis->GetBinLowEdge(ibin);
245  else if (findBinBy == kBinUpEdge)
246  ibinEdge = axis->GetBinUpEdge(ibin);
247  else
248  throw runtime_error("EBinEdge has value " + to_string(findBinBy) + ", which is not defined.");
249  if (abs(ibinEdge - binEdge) < epsilon){
250  foundBin = ibin;
251  break;
252  }
253  }
254  if (foundBin < 1)
255  throw runtime_error("Could not find bin matching " + to_string(binEdge));
256  return foundBin;
257 }
258 
259 // Overlay all generators over data
260 void OverlayGenerators(map<string, map<string, TH1* > > generators, string var, bool drawSpline = false, TLegend ** leg = NULL, double binLowEdge = -1.)
261 {
262  int smoothFactor = 5;
263  cout << "Starting bin checks for var: " << var << endl;
264  for (pair<string, map<string, TH1* > > generator : generators){
265  if (drawSpline)
266  {
267  TSpline3 * spline;
268  if (var == "MuKin"){
269  smoothFactor = 10;
270  int generatorBin = GetBinByEdge(generator.second["MuKin"]->GetXaxis(), binLowEdge, kBinLowEdge);
271  TH2 * hist = (TH2*) generator.second[var];
272  TH1 * histSlice = hist->ProjectionY((hist->GetName() + string("_") + to_string(generatorBin)).c_str(), generatorBin, generatorBin);
273  // cout << "Using bin low edge: " << hist->GetXaxis()->GetBinLowEdge(generatorBin) << " for generator." << endl;
274  histSlice->Rebin(smoothFactor);
275  histSlice->Scale(1.0/double(smoothFactor));
276  // cout << "Generator: " << generator.first << " using generatorBin: " << generatorBin << " (" <<
277  // hist->GetXaxis()->GetBinLowEdge(generatorBin) << " - " << hist->GetXaxis()->GetBinUpEdge(generatorBin) << ")" << endl;
278  spline = new TSpline3(histSlice);
279  }
280  else{
281  generator.second[var]->Rebin(smoothFactor);
282  generator.second[var]->Scale(1.0/smoothFactor);
283  spline = new TSpline3(generator.second[var]);
284  }
285 
286  spline->SetLineColor(generator.second[var]->GetLineColor());
287  spline->SetLineWidth(3);
288  spline->Draw("same");
289  }
290  else
291  {
292  TH1 * hist = generator.second[var];
293  if (var == "MuKin"){
294  int generatorBin = GetBinByEdge(generator.second["MuKin"]->GetXaxis(), binLowEdge, kBinLowEdge);
295  TH2 * hist2D = (TH2*) generator.second[var];
296  hist = hist2D->ProjectionY((hist2D->GetName() + string("_") + to_string(generatorBin)).c_str(), generatorBin, generatorBin);
297  // cout << "Generator: " << generator.first << " using generatorBin: " << generatorBin << " (" <<
298  // hist2D->GetXaxis()->GetBinLowEdge(generatorBin) << " - " << hist2D->GetXaxis()->GetBinUpEdge(generatorBin) << ")" << endl;
299  }
300  hist->SetLineColor(generator.second[var]->GetLineColor());
301  hist->SetLineWidth(3);
302  hist->Draw("hist same");
303  }
304 
305  if (leg)
306  (*leg)->AddEntry(generator.second[var], generator.first.c_str());
307  }
308  if (leg){
309  TLegend * oldLeg = (TLegend*) (*leg)->Clone();
310  (*leg) = ana::AutoPlaceLegend(0.34, 0.2);
311  (*leg)->SetFillStyle(0);
312  (*leg)->SetTextSize(0.03);
313  for (TObject * obj : *(oldLeg->GetListOfPrimitives())){
314  TLegendEntry * entry = (TLegendEntry*) obj;
315  (*leg)->AddEntry(entry->GetObject(), entry->GetLabel());
316  }
317  }
318 }
319 
320 pair<TH1*, TH1*> GeneratorRatio(TH1 * gen, TH1 * data)
321 {
322  TH1 * genClone = (TH1*) gen->Clone();
323  TH1 * genRat = (TH1*) gen->Clone();
324  assert(gen->GetNbinsX() == data->GetNbinsX());
325  genClone->SetTitle(data->GetTitle());
326  genRat->Divide(data);
327  genRat->GetYaxis()->SetTitle("Gen / Data");
328  // genRat->SetLineWidth(2.);
329  return {genClone, genRat};
330 }
331 
332 void PlotGenRatios(map<string, map<string, TH1* > > generators, TH1* dataFull, TH1* dataStat, string var, TLegend ** leg = NULL, double binLowEdge = -1.)
333 {
334  TRatioPlot * baseRatio = NULL;
335  vector<pair<TH1*, TH1*> > ratios;
336  for (pair<string, map<string, TH1* > > generator : generators){
337  TH1 * hist = generator.second[var];
338  if (var == "MuKin"){
339  TH2 * hist2D = (TH2*) generator.second[var];
340  int generatorBin = GetBinByEdge(hist2D->GetXaxis(), binLowEdge, kBinLowEdge);
341  hist = hist2D->ProjectionY((hist->GetName() + string("_") + to_string(generatorBin)).c_str(), generatorBin, generatorBin);
342  hist->SetLineColor(hist2D->GetLineColor());
343  }
344  pair<TH1*, TH1*> genRatio = GeneratorRatio(hist, dataFull);
345  ratios.push_back(genRatio);
346  if (!baseRatio){
347  baseRatio = new TRatioPlot(hist, dataFull);
348  // (*leg)->AddEntry(baseRatio->GetUpperRefObject(), generator.first.c_str());
349  }
350  if(leg)
351  (*leg)->AddEntry(hist, generator.first.c_str());
352  }
353 
354  double maxhist = (*std::max_element(ratios.begin(), ratios.end(), [](auto i, auto j){
355  return (i.first->GetMaximum() < j.first->GetMaximum());
356  })).first->GetMaximum();
357  double maxrat = (*std::max_element(ratios.begin(), ratios.end(), [](auto i, auto j){return i.second->GetMaximum() < j.second->GetMaximum();})).second->GetMaximum();
358  double minrat = (*std::min_element(ratios.begin(), ratios.end(), [](auto i, auto j){return i.second->GetMinimum() < j.second->GetMinimum();})).second->GetMinimum();
359 
360  baseRatio->Draw();
361  baseRatio->SetSplitFraction(0.4);
362  baseRatio->SetSeparationMargin(0.0);
363  baseRatio->GetLowerRefYaxis()->SetTitle("Gen / Data");
364  ((TH1*) baseRatio->GetUpperRefObject())->SetLineWidth(0.0);
365  ((TH1*) baseRatio->GetUpperRefObject())->SetTitle("Generator Ratios");
366  ((TH1*) baseRatio->GetUpperRefObject())->SetTitleSize(0.035, "Y");
367  baseRatio->GetLowerRefGraph()->SetMarkerSize(0.0);
368  baseRatio->GetLowerRefGraph()->GetYaxis()->SetTitleSize(0.035);
369  baseRatio->SetGraphDrawOpt("X");
370  baseRatio->GetUpperRefYaxis()->SetRangeUser(0.0, maxhist * 1.2);
371  baseRatio->GetLowerRefYaxis()->SetRangeUser(minrat * 0.9, 1.1 * maxrat);
372  baseRatio->GetUpperRefXaxis()->SetTitle(dataFull->GetXaxis()->GetTitle());
373  baseRatio->GetLowYaxis()->SetTitle(dataFull->GetXaxis()->GetTitle());
374  baseRatio->GetUpperRefYaxis()->SetTitle(dataFull->GetYaxis()->GetTitle());
375  baseRatio->GetUpperRefYaxis()->CenterTitle(false);
376  baseRatio->GetLowerRefYaxis()->CenterTitle(true);
377  baseRatio->GetLowerRefYaxis()->SetTitleOffset(1.3);
378  baseRatio->SetGridlines({});
379  baseRatio->Draw("noconfint");
380  for(pair<TH1*, TH1*> ratio : ratios)
381  {
382  TH1 * gen = ratio.first;
383  TH1 * rat = ratio.second;
384  baseRatio->GetUpperPad()->cd();
385  gen->Draw("hist same");
386  baseRatio->GetLowerPad()->cd();
387  rat->Draw("hist same");
388  }
389  baseRatio->GetUpperPad()->cd();
390  dataFull->Draw("e1 same");
391  dataStat->Draw("e1 same");
392 
393  if (leg){
394  TLegend * oldLeg = (TLegend*) (*leg)->Clone();
395  (*leg) = ana::AutoPlaceLegend(0.34, var == "MuKin" ? 0.3 : 0.2);
396  (*leg)->SetFillStyle(0);
397  (*leg)->SetTextSize(0.03);
398  for (TObject * obj : *(oldLeg->GetListOfPrimitives())){
399  TLegendEntry * entry = (TLegendEntry*) obj;
400  (*leg)->AddEntry(entry->GetObject(), entry->GetLabel());
401  }
402  (*leg)->Draw();
403  }
404 }
405 
406 // Take a TH2 and print all angle slices
407 void PrintAllSlices(TH2* histSystPlusStat, TH2* histStat, string histLabel, string filename, bool drawSpline, map<string, map<string, TH1* > > generators = {}, bool printRatios = false)
408 {
409  TCanvas * canv = new TCanvas("sliceCanv", "sliceCanv", 1600, 1200);
410  canv->SetLeftMargin(0.14);
411  canv->SetRightMargin(0.05);
412  canv->SetBottomMargin(0.12);
413  TLegend * leg;
414  for (int ibinx = 2; ibinx <= histSystPlusStat->GetNbinsX(); ++ibinx){
415  TH1 * proj0 = histSystPlusStat->ProjectionY((histSystPlusStat->GetName() + string("_") + to_string(ibinx)).c_str(), ibinx, ibinx);
416  proj0->SetTitle(TString::Format("%3.2f < cos(#theta_{#mu}) < %3.2f", histSystPlusStat->GetXaxis()->GetBinLowEdge(ibinx), histSystPlusStat->GetXaxis()->GetBinUpEdge(ibinx)));
417  proj0->GetYaxis()->SetTitle(xsec2DTitle.c_str());
418  proj0->GetXaxis()->SetRangeUser(0.5, 2.5);
419  proj0->GetYaxis()->SetRangeUser(0.0, proj0->GetMaximum() * 1.25);
420  proj0->SetMarkerStyle(histSystPlusStat->GetMarkerStyle());
421  proj0->Draw("E1");
422  leg = new TLegend();
423  leg->SetFillStyle(0);
424  leg->AddEntry(proj0, histLabel.c_str());
425  proj0->Draw("E1 SAME");
426  OverlayGenerators(generators, "MuKin", drawSpline, &leg, histSystPlusStat->GetXaxis()->GetBinLowEdge(ibinx));
427  TH1 * projStat0 = histStat->ProjectionY((histStat->GetName() + string("_") + to_string(ibinx)).c_str(), ibinx, ibinx);
428  projStat0->SetLineColor(kGray);
429  projStat0->Draw("E1 SAME");
430  leg->Draw();
431  canv->SaveAs((filename + "_bin" + to_string(ibinx) + ".pdf").c_str());
432  canv->Clear();
433  if (printRatios)
434  {
435  leg->Clear();
436  leg->AddEntry(proj0, histLabel.c_str());
437  PlotGenRatios(generators, proj0, projStat0, "MuKin", &leg, histSystPlusStat->GetXaxis()->GetBinLowEdge(ibinx));
438  leg->Draw();
439  canv->SaveAs((filename + "_bin" + to_string(ibinx) + "_ratios.pdf").c_str());
440  canv->Clear();
441  }
442  leg->Clear();
443  delete leg;
444  }
445  canv->Clear();
446  delete canv;
447 }
448 
450 {
451  for (int ibinx = 1; ibinx <= hist->GetNbinsX(); ++ibinx){
452  hist->SetBinContent(ibinx, hist->GetBinContent(ibinx) / hist->GetXaxis()->GetBinCenter(ibinx));
453  hist->SetBinError(ibinx, hist->GetBinError(ibinx) / hist->GetXaxis()->GetBinCenter(ibinx));
454  }
455 }
456 
457 map<string, map<string, string> > generator_hist_names
458 {
459  {"genie3",
460  {
461  {"MuKin", "T_muCos_theta"},
462  {"ENu", "xsecVsEnu_2"},
463  {"Q2", "xsecVsq2_2"}
464  }
465  },
466  {"genie3_rebin",
467  {
468  {"MuKin", "T_muCos_theta"},
469  {"ENu", "xsecVsEnu_2"},
470  {"Q2", "xsecVsq2VB_2"}
471  }
472  },
473  {"default",
474  {
475  {"MuKin", "T_muCos_theta_Soup"},
476  {"ENu", "xsecVsEnu_2_Soup"},
477  {"Q2", "xsecVsq2_2_Soup"}
478  }
479  },
480  {"rebin",
481  {
482  {"MuKin", "T_muCos_theta_Soup"},
483  {"ENu", "xsecVsEnu_2_Soup"},
484  {"Q2", "xsecVsq2VB_2_Soup"}
485  }
486  }
487 };
488 
489 // Load the generator histograms from a file, scale if needed
490 map<string, TH1*> LoadGenerator(TFile * file, string label, Color_t col, double scaling = -1, string histLabel = "default")
491 {
492  map<string, TH1*> generatorMap;
493 
494  assert(generator_hist_names.find(histLabel) != generator_hist_names.end());
495  assert(generator_hist_names[histLabel].find("MuKin") != generator_hist_names[histLabel].end());
496  assert(generator_hist_names[histLabel].find("ENu") != generator_hist_names[histLabel].end());
497  assert(generator_hist_names[histLabel].find("Q2") != generator_hist_names[histLabel].end());
498 
499  generatorMap["MuKin"] = (TH2*) file->Get(generator_hist_names[histLabel]["MuKin"].c_str())->Clone();
500  generatorMap["ENu"] = (TH1*) file->Get(generator_hist_names[histLabel]["ENu"].c_str())->Clone();
501  generatorMap["Q2"] = (TH1*) file->Get(generator_hist_names[histLabel]["Q2"].c_str())->Clone();
502 
503  generatorMap["MuKin"]->SetLineColor(col);
504  generatorMap["ENu"] ->SetLineColor(col);
505  generatorMap["Q2"] ->SetLineColor(col);
506 
507  string appendLabel = histLabel == "rebin" || histLabel == "genie3_rebin" ? "rebin" : "";
508  generatorMap["MuKin"]->SetName((label + "_MuKin" + appendLabel).c_str());
509  generatorMap["ENu"] ->SetName((label + "_ENu" + appendLabel).c_str());
510  generatorMap["Q2"] ->SetName((label + "_Q2" + appendLabel).c_str());
511 
512  if (scaling > 0)
513  {
514  generatorMap["MuKin"]->Scale(scaling);
515  generatorMap["ENu"] ->Scale(scaling);
516  generatorMap["Q2"] ->Scale(scaling);
517  }
518 
519  generatorMap["MuKin"]->GetXaxis()->SetRangeUser(0.5, 1.0);
520  generatorMap["MuKin"]->GetYaxis()->SetRangeUser(0.5, 2.5);
521  generatorMap["ENu"] ->GetXaxis()->SetRangeUser(0.5, 5.0);
522  generatorMap["Q2"] ->GetXaxis()->SetRangeUser(0.0, 2.8);
523 
524  // Set XSec of ENu to XSec / ENu of ENu
525  SigmaDivideENu(generatorMap["ENu"]);
526 
527  return generatorMap;
528 }
529 
530 map<string, TH1*> LoadGeniePrediction(TFile * file, string label, Color_t col, string weight)
531 {
532  map<string, TH1*> geniePred;
533  geniePred["MuKin"] = (TH2*) file->Get(("geniePrediction/" + weight + "/MuKin").c_str())->Clone();
534  geniePred["ENu"] = (TH1*) file->Get(("geniePrediction/" + weight + "/ENu").c_str())->Clone();
535  geniePred["Q2"] = (TH1*) file->Get(("geniePrediction/" + weight + "/Q2").c_str())->Clone();
536  geniePred["MuKin"]->SetLineColor(col);
537  geniePred["ENu"] ->SetLineColor(col);
538  geniePred["Q2"] ->SetLineColor(col);
539  geniePred["MuKin"]->SetName((label + "_MuKin").c_str());
540  geniePred["ENu"] ->SetName((label + "_ENu" ).c_str());
541  geniePred["Q2"] ->SetName((label + "_Q2" ).c_str());
542 
543  // Get rid of huge scaling factor
544  geniePred["MuKin"]->Scale(1E39);
545  geniePred["ENu"] ->Scale(1E39);
546  geniePred["Q2"] ->Scale(1E39);
547 
548  geniePred["MuKin"]->GetXaxis()->SetRangeUser(0.5, 1.0);
549  geniePred["MuKin"]->GetYaxis()->SetRangeUser(0.5, 2.5);
550  geniePred["ENu"] ->GetXaxis()->SetRangeUser(0.5, 5.0);
551  geniePred["Q2"] ->GetXaxis()->SetRangeUser(0.0, 2.8);
552 
553  // Set XSec of ENu to XSec / ENu of ENu
554  SigmaDivideENu(geniePred["ENu"]);
555 
556  return geniePred;
557 }
558 
559 void AreaNormalize(TH1* generator, TH1* data, const vector<double>& binning, const string label)
560 {
561  double genBinLow = generator->FindBin(binning[0]);
562  double genBinUp = generator->FindBin(binning[binning.size() - 1]) - 1;
563  double dataBinLow = data->FindBin(binning[0]);
564  double dataBinUp = data->FindBin(binning[binning.size() - 1]) - 1;
565  cout << "Scaling " << label << ": " << data->Integral(dataBinLow, dataBinUp, "width") / generator->Integral(genBinLow, genBinUp, "width") << endl;
566  generator->Scale(data->Integral(dataBinLow, dataBinUp, "width") / generator->Integral(genBinLow, genBinUp, "width"));
567 }
568 
569 // TH2 is slightly different. We need to normalize based on only the bins we present, not outside the phase space
570 void AreaNormalize2D(TH2* generator, TH2* data, const string label)
571 {
572  double genIntegral = 0.;
573  double dataIntegral = 0.;
574  for(int ibinx = 1; ibinx <= data->GetNbinsX(); ++ibinx)
575  for(int ibiny = 1; ibiny <= data->GetNbinsY(); ++ibiny)
576  if (data->GetBinContent(ibinx, ibiny) > 0.0){
577  // cout << "\tFound (" << std::scientific << to_string(ibinx) << "," << to_string(ibiny) << ") in data.\n\tAdding: "
578  // << data->GetBinContent(ibinx, ibiny) << " * " << to_string(data->GetXaxis()->GetBinWidth(ibinx)) << " * " << to_string(data->GetYaxis()->GetBinWidth(ibiny)) << endl;
579  dataIntegral += data->GetBinContent(ibinx, ibiny) * data->GetXaxis()->GetBinWidth(ibinx) * data->GetYaxis()->GetBinWidth(ibiny);
580  }
581  for(int ibinx = 1; ibinx <= generator->GetNbinsX(); ++ibinx)
582  for(int ibiny = 1; ibiny <= generator->GetNbinsY(); ++ibiny){
583  int dataBinX = data->GetXaxis()->FindBin(generator->GetXaxis()->GetBinCenter(ibinx));
584  int dataBinY = data->GetYaxis()->FindBin(generator->GetYaxis()->GetBinCenter(ibiny));
585  if (data->GetBinContent(dataBinX, dataBinY) > 0.0){
586  // cout << "\tFound (" << std::scientific << to_string(ibinx) << "," << to_string(ibiny) << ") in generator.\n\tAdding: "
587  // << generator->GetBinContent(ibinx, ibiny) << " * " << to_string(generator->GetXaxis()->GetBinWidth(ibinx)) << " * " << to_string(generator->GetYaxis()->GetBinWidth(ibiny)) << endl;
588  genIntegral += generator->GetBinContent(ibinx, ibiny) * generator->GetXaxis()->GetBinWidth(ibinx) * generator->GetYaxis()->GetBinWidth(ibiny);
589  }
590  }
591  cout << "Scaling: " << label << ": " << dataIntegral / genIntegral << endl;
592  generator->Scale(dataIntegral / genIntegral);
593 }
594 
595 int Get2DBinArea(TH2* hist, int binxLow, int binxHigh, int binyLow, int binyHigh)
596 {
597  double widthx = 0;
598  double widthy = 0;
599  for (int ibinx = binxLow; ibinx <= binxHigh; ++ibinx)
600  widthx += hist->GetXaxis()->GetBinWidth(ibinx);
601  for (int ibiny = binyLow; ibiny <= binyHigh; ++ibiny)
602  widthy += hist->GetYaxis()->GetBinWidth(ibiny);
603  return widthx * widthy;
604 }
605 
606 map<string, map<string, TH1* > > RebinGenerators(map<string, map<string, TH1* > > generators)
607 {
608  map<string, map<string, TH1* > > rebinnedGenerators;
609  for (pair<string, map<string, TH1* > > generatorPair : generators){
610  string genLabel = generatorPair.first;
611  map<string, TH1* > & generator = generatorPair.second;
612 
613  // Initialize the empty histograms to the right binning
614  TH2D * mukin_hist = new TH2D((generator["MuKin"]->GetName() + string("_databins")).c_str(), (genLabel + " Rebin").c_str(),
616  TH1D * enu_hist = new TH1D((generator["ENu"]->GetName() + string("_databins")).c_str(), (genLabel + " Rebin").c_str(),
617  enuBinEdges.size() - 1, enuBinEdges.data());
618  TH1D * q2_hist = new TH1D((generator["Q2"]->GetName() + string("_databins")).c_str(), (genLabel + " Rebin").c_str(),
619  q2BinEdges.size() - 1, q2BinEdges.data());
620 
621  mukin_hist->SetLineColor(generator["MuKin"]->GetLineColor());
622  enu_hist ->SetLineColor(generator["ENu"] ->GetLineColor());
623  q2_hist ->SetLineColor(generator["Q2"]->GetLineColor());
624 
625  // Loop through MuKin bin edges, fill with correct integral, scale by inverse bin width.
626  for (unsigned int ibinx = 0; ibinx < angBinEdges_reduced.size() - 1; ++ibinx)
627  {
628  double binxEdgeLow = angBinEdges_reduced[ibinx];
629  double binxEdgeUp = angBinEdges_reduced[ibinx + 1];
630  int binxLow = GetBinByEdge(generator["MuKin"]->GetXaxis(), binxEdgeLow, kBinLowEdge);
631  int binxHigh = GetBinByEdge(generator["MuKin"]->GetXaxis(), binxEdgeUp, kBinUpEdge);
632  for (unsigned int ibiny = 0; ibiny < mukeBinEdges_reduced.size() - 1; ++ibiny)
633  {
634  double binyEdgeLow = mukeBinEdges_reduced[ibiny];
635  double binyEdgeUp = mukeBinEdges_reduced[ibiny + 1];
636  int binyLow = GetBinByEdge(generator["MuKin"]->GetYaxis(), binyEdgeLow, kBinLowEdge);
637  int binyHigh = GetBinByEdge(generator["MuKin"]->GetYaxis(), binyEdgeUp, kBinUpEdge);
638  double binIntegral = ((TH2*) generator["MuKin"])->Integral(binxLow, binxHigh, binyLow, binyHigh);
639  int numBins = (binxHigh - binxLow + 1) * (binyHigh - binyLow + 1);
640  mukin_hist->SetBinContent(ibinx + 1, ibiny + 1, binIntegral / numBins);
641  }
642  }
643 
644  // Simpler logic now for the ENu
645  for (unsigned int ibinx = 0; ibinx < enuBinEdges.size() - 1; ++ibinx)
646  {
647  double binEdgeLow = enuBinEdges[ibinx];
648  double binEdgeUp = enuBinEdges[ibinx + 1];
649  int binLow = GetBinByEdge(generator["ENu"]->GetXaxis(), binEdgeLow, kBinLowEdge);
650  int binHigh = GetBinByEdge(generator["ENu"]->GetXaxis(), binEdgeUp, kBinUpEdge);
651  double binIntegral = generator["ENu"]->Integral(binLow, binHigh);
652  enu_hist->SetBinContent(ibinx + 1, binIntegral / (binHigh - binLow + 1));
653  }
654 
655  // Fill the Q2
656  for (unsigned int ibinx = 0; ibinx < q2BinEdges.size() - 1; ++ibinx)
657  {
658  double binEdgeLow = q2BinEdges[ibinx];
659  double binEdgeUp = q2BinEdges[ibinx + 1];
660  int binLow = GetBinByEdge(generator["Q2"]->GetXaxis(), binEdgeLow, kBinLowEdge);
661  int binHigh = GetBinByEdge(generator["Q2"]->GetXaxis(), binEdgeUp, kBinUpEdge);
662  double binIntegral = generator["Q2"]->Integral(binLow, binHigh);
663  q2_hist->SetBinContent(ibinx + 1, binIntegral / (binHigh - binLow + 1));
664  }
665 
666  rebinnedGenerators[generatorPair.first]["MuKin"] = mukin_hist;
667  rebinnedGenerators[generatorPair.first]["ENu"] = enu_hist;
668  rebinnedGenerators[generatorPair.first]["Q2"] = q2_hist;
669  }
670  return rebinnedGenerators;
671 }
672 
673 // Chi2 and minimization
674 
675 double CalculateENuChi2(const TMatrixD& invCovMat, TH1 * data, TH1 * generator, string label)
676 {
677  double chi2 = 0;
678  double naiveChi2 = 0;
679  assert(invCovMat.GetNrows() == invCovMat.GetNcols());
680  assert(invCovMat.GetNrows() == (int) enuBinEdges.size() - 1);
681  assert(invCovMat.GetNrows() == data->GetNbinsX());
682  assert(invCovMat.GetNrows() == generator->GetNbinsX());
683  TH2D * chi2Contribs = new TH2D((label + "_enu_chi2_contribs").c_str(), (label + " E_{#nu} #chi^{2} Contributions").c_str(), invCovMat.GetNrows(), 0, invCovMat.GetNrows(), invCovMat.GetNcols(), 0, invCovMat.GetNcols());
684  for (int ibin = 0; ibin < invCovMat.GetNrows(); ++ibin)
685  for (int jbin = 0; jbin < invCovMat.GetNcols(); ++jbin)
686  {
687  double binChi2 = (generator->GetBinContent(ibin + 1) - data->GetBinContent(ibin + 1)) * invCovMat[ibin][jbin] * (generator->GetBinContent(jbin + 1) - data->GetBinContent(jbin + 1));
688  chi2Contribs->SetBinContent(ibin + 1, jbin + 1, binChi2);
689  chi2 += binChi2;
690  if (ibin == jbin) // Only diagonals for naive calculation
691  naiveChi2 += pow(generator->GetBinContent(ibin + 1) - data->GetBinContent(ibin + 1), 2) / pow(data->GetBinError(ibin + 1), 2);
692  }
693  chi2Contribs->Draw("COLZ");
694  cout << "\tNaive Chi2: " << naiveChi2 << endl;
695  return chi2;
696 }
697 
698 double Calculate1DChi2_simple(const TMatrixD& invCovMat, const TH1 * data, const TH1 * generator)
699 {
700  double chi2 = 0;
701  for (int ibin = 0; ibin < invCovMat.GetNrows(); ++ibin)
702  for (int jbin = 0; jbin < invCovMat.GetNcols(); ++jbin)
703  chi2 += (generator->GetBinContent(ibin + 1) - data->GetBinContent(ibin + 1)) * invCovMat[ibin][jbin] * (generator->GetBinContent(jbin + 1) - data->GetBinContent(jbin + 1));
704  return chi2;
705 }
706 
710 void fcn_enu(int &npar, double *gin, double &f, double *par, int iflag)
711 {
712  cout << "Running fcn_enu() at scaling: " << to_string(par[0]) << endl;
713  TH1 * scaled_gen = (TH1*) chi2_min_enu_gen->Clone("temp");
714  scaled_gen->Scale(par[0]);
715  f = Calculate1DChi2_simple(*chi2_min_enu_invCovMat, chi2_min_enu_data, scaled_gen);
716  cout << "\tGives chi2: " << to_string(f) << endl;
717  scaled_gen->Delete();
718 }
719 
720 double CalculateQ2Chi2(const TMatrixD& invCovMat, TH1 * data, TH1 * generator, string label)
721 {
722  double chi2 = 0;
723  double naiveChi2 = 0;
724  assert(invCovMat.GetNrows() == invCovMat.GetNcols());
725  assert(invCovMat.GetNrows() == (int) q2BinEdges.size() - 1);
726  assert(invCovMat.GetNrows() == data->GetNbinsX());
727  assert(invCovMat.GetNrows() == generator->GetNbinsX());
728  TH2D * chi2Contribs = new TH2D((label + "_q2_chi2_contribs").c_str(), (label + " q^{2} #chi^{2} Contributions").c_str(), invCovMat.GetNrows(), 0, invCovMat.GetNrows(), invCovMat.GetNcols(), 0, invCovMat.GetNcols());
729  for (int ibin = 0; ibin < invCovMat.GetNrows(); ++ibin)
730  for (int jbin = 0; jbin < invCovMat.GetNcols(); ++jbin)
731  {
732  double binChi2 = (generator->GetBinContent(ibin + 1) - data->GetBinContent(ibin + 1)) * invCovMat[ibin][jbin] * (generator->GetBinContent(jbin + 1) - data->GetBinContent(jbin + 1));
733  chi2Contribs->SetBinContent(ibin + 1, jbin + 1, binChi2);
734  chi2 += binChi2;
735  if (ibin == jbin) // Only diagonals for naive calculation
736  naiveChi2 += pow(generator->GetBinContent(ibin + 1) - data->GetBinContent(ibin + 1), 2) / pow(data->GetBinError(ibin + 1), 2);
737  }
738  chi2Contribs->Draw("COLZ");
739  cout << "\tNaive Chi2: " << naiveChi2 << endl;
740  return chi2;
741 }
742 
746 void fcn_q2(int &npar, double *gin, double &f, double *par, int iflag)
747 {
748  cout << "Running fcn_q2() at scaling: " << to_string(par[0]) << endl;
749  TH1 * scaled_gen = (TH1*) chi2_min_q2_gen->Clone("temp");
750  scaled_gen->Scale(par[0]);
751  f = Calculate1DChi2_simple(*chi2_min_q2_invCovMat, chi2_min_q2_data, scaled_gen);
752  cout << "\tGives chi2: " << to_string(f) << endl;
753  scaled_gen->Delete();
754 }
755 
756 double Calculate2DChi2(const TMatrixD& invCovMat, TH2 * data, TH2 * generator, TH1 * covToGlobalBinMap, string label)
757 {
758  double chi2 = 0;
759  double naiveChi2 = 0;
760  assert(invCovMat.GetNrows() == invCovMat.GetNcols());
761  assert(invCovMat.GetNrows() == covToGlobalBinMap->GetNbinsX());
762  assert(data->GetNbinsX() == generator->GetNbinsX());
763  assert(data->GetNbinsY() == generator->GetNbinsY());
764  TH2D * chi2Contribs = new TH2D((label + "_mukin_chi2_contribs").c_str(), (label + "Muon Kinematic #chi^{2} Contributions").c_str(), invCovMat.GetNrows(), 0, invCovMat.GetNrows(), invCovMat.GetNcols(), 0, invCovMat.GetNcols());
765  // Check the binnings are exact with edges
766  for (int ibinx = 1; ibinx <= generator->GetNbinsX(); ++ibinx){
767  for (int ibiny = 1; ibiny <= generator->GetNbinsY(); ++ibiny){
768  double xval = generator->GetXaxis()->GetBinLowEdge(ibinx);
769  double yval = generator->GetYaxis()->GetBinLowEdge(ibiny);
770  if (find_if(angBinEdges.begin(), angBinEdges.end(), [xval](double b) {return abs(xval - b) < 1E-10;}) == angBinEdges.end() ||
771  find_if(mukeBinEdges.begin(), mukeBinEdges.end(), [yval](double b) {return abs(yval - b) < 1E-10;}) == mukeBinEdges.end()){
772  cerr << "\tDid not find bin: (" << generator->GetXaxis()->GetBinLowEdge(ibinx) << ", " << generator->GetYaxis()->GetBinLowEdge(ibiny) << ") in list of bins." << endl;
773  cerr << "\tBinning (2D) does not match! Skipping generator." << endl;
774  return -1;
775  }
776  }
777  }
778  for (int ibin = 0; ibin < invCovMat.GetNrows(); ++ibin)
779  for (int jbin = 0; jbin < invCovMat.GetNcols(); ++jbin)
780  {
781  int iGlobalBin = covToGlobalBinMap->GetBinContent(ibin + 1);
782  int jGlobalBin = covToGlobalBinMap->GetBinContent(jbin + 1);
783  int ibinX = -1;
784  int ibinY = -1;
785  int jbinX = -1;
786  int jbinY = -1;
787  int temp = -1;
788  data->GetBinXYZ(iGlobalBin, ibinX, ibinY, temp); // Fill the binx and biny for global bin I
789  data->GetBinXYZ(jGlobalBin, jbinX, jbinY, temp); // Fill the binx and biny for global bin J
790  double binChi2 = (generator->GetBinContent(ibinX, ibinY) - data->GetBinContent(ibinX, ibinY)) * invCovMat[ibin][jbin] * (generator->GetBinContent(jbinX, jbinY) - data->GetBinContent(jbinX, jbinY));
791  chi2Contribs->SetBinContent(ibin + 1, jbin + 1, binChi2);
792  if (ibin == jbin) // Only diagonals for naive calculation
793  naiveChi2 += pow(generator->GetBinContent(ibinX, ibinY) - data->GetBinContent(ibinX, ibinY), 2) / pow(data->GetBinError(ibinX, ibinY), 2);
794  // cout << "\tUsing matrix bins (" << ibin << "," << jbin << ")" << endl;
795  // cout << "\tCorresponding to MuKin:((" << ibinXCenter << ", " << ibinYCenter << "), (" << jbinXCenter << ", " << jbinYCenter << "))" << endl;
796  // cout << "\tUsing generator bins ((" << ibinXGen << "," << ibinYGen << "), (" << jbinXGen << ", " << jbinYGen << "))" << endl;
797  // cout << "\tUsing data bins ((" << ibinXData << "," << ibinYData << "), (" << jbinXData << ", " << jbinYData << "))" << endl;
798  // cout << "\tAdding (" << generator->GetBinContent(ibinXGen, ibinYGen) << " - " << data->GetBinContent(ibinXData, ibinYData) << ") * " <<
799  // invCovMat[ibin][jbin] << " * (" << generator->GetBinContent(jbinXGen, jbinYGen) << " - " << data->GetBinContent(jbinXData, jbinYData) << ")" << endl;
800  // cout << "\t(=" << binChi2 << ")" << endl;
801  chi2 += binChi2;
802  }
803  chi2Contribs->Draw("COLZ");
804  cout << "\tNaive Chi2: " << naiveChi2 << endl;
805  return chi2;
806 }
807 
808 void PlotDataGenerators(TH1 * dataSystPlusStat, TH1 * dataStat, string dataTitle, map<string, map<string, TH1* > > &generatorMap, string varName, bool drawSpline = false)
809 {
810  dataSystPlusStat->Draw("E1");
811  TLegend * leg = new TLegend(0.6, 0.6, 0.9, 0.85);
812  leg->SetFillStyle(0);
813  leg->AddEntry(dataSystPlusStat, dataTitle.c_str());
814  OverlayGenerators(generatorMap, varName, drawSpline, &leg);
815  dataSystPlusStat->Draw("E1 SAME");
816  dataStat->SetLineColor(kGray);
817  dataStat->Draw("E1 SAME");
818  leg->Draw();
819 }
820 
822 {
823  int nrowcols = hist->GetNbinsX();
824  TMatrixD statErrs(nrowcols, nrowcols);
825  for(int irow = 0; irow < nrowcols; irow++)
826  statErrs[irow][irow] = hist->GetBinError(irow+1);
827  return statErrs;
828 }
829 
830 void PlotCovMatrices(TMatrixD & cov, TMatrixD & inv, string outputFolder, string label, string title)
831 {
832  TCanvas * tempCanv = new TCanvas("temp", "temp", 1600, 1200);
833  tempCanv->SetRightMargin(0.15);
834 
835  TH2D * covHist = new TH2D(cov);
836  covHist->SetName((label + "CovHist").c_str());
837  covHist->SetTitle((title + " Covariance Matrix;;;#times 10^{-78}").c_str());
838  covHist->Draw("COLZ");
839  tempCanv->SaveAs((outputFolder + "/" + label + "CovMat.pdf").c_str());
840 
841  TH2D * invHist = new TH2D(inv);
842  invHist->SetName((label + "InvCovHist").c_str());
843  invHist->SetTitle((title + " Inv. Cov Matrix;;;#times 10^{78}").c_str());
844  invHist->Draw("COLZ");
845  tempCanv->SaveAs((outputFolder + "/" + label + "InvCovMat.pdf").c_str());
846 
847  TMatrixD identity = cov * inv;
848  TH2D * identityHist = new TH2D(identity);
849  identityHist->SetName((label + "Identity").c_str());
850  identityHist->SetTitle((title + " Cov #times InvCov").c_str());
851  identityHist->Draw("COLZ");
852  tempCanv->SaveAs((outputFolder + "/" + label + "Identity.pdf").c_str());
853 
854  covHist->Write();
855  invHist->Write();
856 
857  tempCanv->Clear();
858  delete tempCanv;
859 }
860 
861 void PlotCovariance(string outputFolder = ".")
862 {
863  gStyle->SetTitleYSize(0.05);
864  // gStyle->SetTitleZSize(0.15);
865 
866  TCanvas * c1 = new TCanvas("c1", "c1", 1600, 1200);
867  string shapeOnlyFilename = "CovCurrent/Calculate_covariances_shape_all.root";
868  string totalErrFilename = "CovCurrent/Calculate_covariances_total_all.root";
869 
870  TFile * outFile = TFile::Open((outputFolder + "/numuCCXSecResults.root").c_str(), "RECREATE");
871  TDirectory * xsecDir = outFile->mkdir("results");
872  TDirectory * generatorDir = outFile->mkdir("generators");
873 
874  // Get a 2D bin map of x-y to
875  TH2I * bin2DToCov = Make2DToCovBin();
876 
877  /////////////////////////////////
878  // Load Statistical Errors
879  /////////////////////////////////
880  TFile * Stat2DFile = TFile::Open("/nova/ana/users/connorj/NumuCC/NERSC_Run_Dec2019/xsecUncertResults.root");
881  TFile * Stat1DFile = TFile::Open("/nova/ana/users/connorj/NumuCC/Feb2020_MuonE_Angle_SystFix/xsecUncertResults_1DOnly.root");
882  assert(Stat2DFile && Stat2DFile->IsOpen());
883  assert(Stat1DFile && Stat1DFile->IsOpen());
884  TH2 * statMuKin_originalbins = (TH2*) Stat2DFile->Get("cv/xsec_cv_2DProj");
885  TH1 * statENu_originalbins = (TH1*) Stat1DFile->Get("cv/xsec_cv_ENu");
886  TH1 * statQ2_originalbins = (TH1*) Stat1DFile->Get("cv/xsec_cv_Q2");
887  statMuKin_originalbins->Scale(1E39);
888  statENu_originalbins->Scale(1E39);
889  statQ2_originalbins->Scale(1E39);
890 
891  TH2 * statMuKin = new TH2D("statMuKin", ("Statistical Errors;" + muAngTitle + ";" + mukeTitle + ";" + xsec2DTitle).c_str(),
892  angBinEdges_reduced.size() - 1, angBinEdges_reduced.data(),
893  mukeBinEdges_reduced.size() - 1, mukeBinEdges_reduced.data());
894  TH1 * statENu = statENu_originalbins->Rebin(enuBinEdges.size() - 1, "statENu", enuBinEdges.data());
895  TH1 * statQ2 = statQ2_originalbins ->Rebin(q2BinEdges.size() - 1, "statQ2", q2BinEdges.data());
896  // Loop through MuKin bin edges, fill with correct integral, scale by inverse bin width.
897  for (unsigned int ibinx = 0; ibinx < angBinEdges_reduced.size() - 1; ++ibinx)
898  {
899  double binxEdgeLow = angBinEdges_reduced[ibinx];
900  double binxEdgeUp = angBinEdges_reduced[ibinx + 1];
901  int binxLow = GetBinByEdge(statMuKin_originalbins->GetXaxis(), binxEdgeLow, kBinLowEdge);
902  int binxHigh = GetBinByEdge(statMuKin_originalbins->GetXaxis(), binxEdgeUp, kBinUpEdge);
903  for (unsigned int ibiny = 0; ibiny < mukeBinEdges_reduced.size() - 1; ++ibiny)
904  {
905  double binyEdgeLow = mukeBinEdges_reduced[ibiny];
906  double binyEdgeUp = mukeBinEdges_reduced[ibiny + 1];
907  int binyLow = GetBinByEdge(statMuKin_originalbins->GetYaxis(), binyEdgeLow, kBinLowEdge);
908  int binyHigh = GetBinByEdge(statMuKin_originalbins->GetYaxis(), binyEdgeUp, kBinUpEdge);
909  statMuKin->SetBinContent(ibinx + 1, ibiny + 1, statMuKin_originalbins->GetBinContent(binxLow, binyLow));
910  statMuKin->SetBinError(ibinx + 1, ibiny + 1, statMuKin_originalbins->GetBinError(binxLow, binyLow));
911  }
912  }
913 
914  // Add the statistical errors to the diagonals of covariance matrix
915  TDirectory * statErrDir = xsecDir->mkdir("statErrs");
916  statErrDir->cd();
917  statMuKin->Write();
918  statENu->Write();
919  statQ2->Write();
920 
921  TMatrixD statMatMuKin(bin2DToCov->GetMaximum(), bin2DToCov->GetMaximum());
922  for(int ibin = 0; ibin < statMatMuKin.GetNrows(); ++ibin){
923  unsigned int binx = 0;
924  unsigned int biny = 0;
925  for(int ibinx = 1; ibinx <= bin2DToCov->GetNbinsX(); ++ibinx)
926  for(int ibiny = 1; ibiny <= bin2DToCov->GetNbinsY(); ++ibiny)
927  if (ibin == bin2DToCov->GetBinContent(ibinx, ibiny))
928  {
929  binx = ibinx;
930  biny = ibiny;
931  }
932  assert(binx > 0);
933  assert(biny > 0);
934  statMatMuKin[ibin][ibin] = statMuKin->GetBinError(binx, biny);
935  }
936  TMatrixD statMatENu(StatMatFromHist(statENu));
937  TMatrixD statMatQ2 (StatMatFromHist(statQ2));
938 
939  /////////////////////////////////
940  // Load Total Error
941  /////////////////////////////////
942 
943  // Load the nominal hists
944  TFile * totalErrFile = TFile::Open(totalErrFilename.c_str());
945  TH1 * nomMuKinTotal1D = (TH1*) totalErrFile->Get("hcv_mukin");
946  TH1 * nomENuTotal_nobininfo = (TH1*) totalErrFile->Get("hcv_enu");
947  TH1 * nomQ2Total_nobininfo = (TH1*) totalErrFile->Get("hcv_q2");
948  nomMuKinTotal1D ->Scale(1E39);
949  nomENuTotal_nobininfo->Scale(1E39);
950  nomQ2Total_nobininfo ->Scale(1E39);
951 
952  // Load the covariance matrices
953  TH2 * covTotalMuKin = (TH2*) totalErrFile->Get("cov_mukin");
954  TH2 * covTotalENu = (TH2*) totalErrFile->Get("cov_enu");
955  TH2 * covTotalQ2 = (TH2*) totalErrFile->Get("cov_q2");
956  covTotalMuKin->Scale(1E78);
957  covTotalENu ->Scale(1E78);
958  covTotalQ2 ->Scale(1E78);
959 
960  // Make matricies and inverses
961  TMatrixD covMatTotalMuKin(covTotalMuKin->GetNbinsX(), covTotalMuKin->GetNbinsY());
962  TMatrixD covMatTotalENu (covTotalENu ->GetNbinsX(), covTotalENu ->GetNbinsY());
963  TMatrixD covMatTotalQ2 (covTotalQ2 ->GetNbinsX(), covTotalQ2 ->GetNbinsY());
964  for (int ibinx = 1; ibinx <= covTotalMuKin->GetNbinsX(); ++ibinx)
965  for (int ibiny = 1; ibiny <= covTotalMuKin->GetNbinsY(); ++ibiny)
966  covMatTotalMuKin[ibinx - 1][ibiny - 1] = covTotalMuKin->GetBinContent(ibinx, ibiny);
967  for (int ibinx = 1; ibinx <= covTotalENu->GetNbinsX(); ++ibinx)
968  for (int ibiny = 1; ibiny <= covTotalENu->GetNbinsY(); ++ibiny)
969  covMatTotalENu[ibinx - 1][ibiny - 1] = covTotalENu->GetBinContent(ibinx, ibiny);
970  for (int ibinx = 1; ibinx <= covTotalQ2->GetNbinsX(); ++ibinx)
971  for (int ibiny = 1; ibiny <= covTotalQ2->GetNbinsY(); ++ibiny)
972  covMatTotalQ2[ibinx - 1][ibiny - 1] = covTotalQ2->GetBinContent(ibinx, ibiny);
973  TMatrixD covMatTotalMuKinFull = covMatTotalMuKin + statMatMuKin;
974  TMatrixD covMatTotalENuFull = covMatTotalENu + statMatENu;
975  TMatrixD covMatTotalQ2Full = covMatTotalQ2 + statMatQ2;
976  TMatrixD invCovMatTotalMuKin = covMatTotalMuKinFull;
977  TMatrixD invCovMatTotalENu = covMatTotalENuFull;
978  TMatrixD invCovMatTotalQ2 = covMatTotalQ2Full;
979  invCovMatTotalMuKin.Invert();
980  invCovMatTotalENu .Invert();
981  invCovMatTotalQ2 .Invert();
982 
983  // Set nominal error bars to the sqrt of covariance diagonals
984  SetErrorFromCovDiag(nomMuKinTotal1D , covMatTotalMuKinFull);
985  SetErrorFromCovDiag(nomENuTotal_nobininfo, covMatTotalENuFull);
986  SetErrorFromCovDiag(nomQ2Total_nobininfo , covMatTotalQ2Full);
987 
988  // Load the right binnings from the nominal
989  TH2 * nomMuKinTotal = Construct2D(nomMuKinTotal1D);
990  TH1 * nomENuTotal = ConstructENu(nomENuTotal_nobininfo);
991  TH1 * nomQ2Total = ConstructQ2(nomQ2Total_nobininfo);
992  SigmaDivideENu(nomENuTotal); // ENu then gets scaled.
993 
994  // Load a second copy with the stats errors included
995  TH1 * nomMuKinTotalStats1D = (TH1*) nomMuKinTotal1D ->Clone("hmukin_cv_2D_stats");
996  TH1 * nomENuTotalStats_nobininfo = (TH1*) nomENuTotal_nobininfo->Clone("henu_cv_enuaxis_stats");
997  TH1 * nomQ2TotalStats_nobininfo = (TH1*) nomQ2Total_nobininfo ->Clone("hq2_cv_q2axis_stats");
998  SetErrorFromCovDiag(nomMuKinTotalStats1D , statMatMuKin);
999  SetErrorFromCovDiag(nomENuTotalStats_nobininfo, statMatENu);
1000  SetErrorFromCovDiag(nomQ2TotalStats_nobininfo , statMatQ2);
1001  TH2 * nomMuKinStats = Construct2D(nomMuKinTotalStats1D);
1002  TH1 * nomENuStats = ConstructENu(nomENuTotalStats_nobininfo);
1003  TH1 * nomQ2Stats = ConstructQ2(nomQ2TotalStats_nobininfo);
1004  SigmaDivideENu(nomENuStats); // ENu then gets scaled.
1005 
1006  // Set titles, axis, markers, etc.
1007  nomMuKinTotal->SetTitle((";" + muAngTitle + ";" + mukeTitle + ";" + xsec2DTitle).c_str());
1008  nomENuTotal-> SetTitle((";" + enuTitle + ";" + xsecENuTitle).c_str());
1009  nomQ2Total-> SetTitle((";" + q2Title + ";" + xsecQ2Title).c_str());
1010  nomMuKinTotal->SetMarkerStyle(20);
1011  nomENuTotal-> SetMarkerStyle(20);
1012  nomQ2Total-> SetMarkerStyle(20);
1013 
1014  // Save to the output
1015  TDirectory * totalUncertResults = xsecDir->mkdir("totalUncertResults");
1016  totalUncertResults->cd();
1017  nomMuKinTotal->Write("MuKin");
1018  nomENuTotal->Write("ENu");
1019  nomQ2Total->Write("Q2");
1020 
1021  /////////////////////////////////
1022  // Load Shape-only Error
1023  /////////////////////////////////
1024 
1025  // Load the nominal hists
1026  TFile * shapeOnlyFile = TFile::Open(shapeOnlyFilename.c_str());
1027  TH1 * nomMuKinShape1D = (TH1*) shapeOnlyFile->Get("hcv_mukin");
1028  TH1 * nomENuShape_nobininfo = (TH1*) shapeOnlyFile->Get("hcv_enu");
1029  TH1 * nomQ2Shape_nobininfo = (TH1*) shapeOnlyFile->Get("hcv_q2");
1030  nomMuKinShape1D ->Scale(1E39);
1031  nomENuShape_nobininfo->Scale(1E39);
1032  nomQ2Shape_nobininfo ->Scale(1E39);
1033 
1034  // Load the covariance matrices
1035  TH2 * covShapeMuKin = (TH2*) shapeOnlyFile->Get("cov_mukin");
1036  TH2 * covShapeENu = (TH2*) shapeOnlyFile->Get("cov_enu");
1037  TH2 * covShapeQ2 = (TH2*) shapeOnlyFile->Get("cov_q2");
1038  covShapeMuKin->Scale(1E78);
1039  covShapeENu ->Scale(1E78);
1040  covShapeQ2 ->Scale(1E78);
1041 
1042  // Make matricies and inverseShape
1043  TMatrixD covMatShapeMuKin(covShapeMuKin->GetNbinsX(), covShapeMuKin->GetNbinsY());
1044  TMatrixD covMatShapeENu (covShapeENu ->GetNbinsX(), covShapeENu ->GetNbinsY());
1045  TMatrixD covMatShapeQ2 (covShapeQ2 ->GetNbinsX(), covShapeQ2 ->GetNbinsY());
1046  for (int ibinx = 1; ibinx <= covShapeMuKin->GetNbinsX(); ++ibinx)
1047  for (int ibiny = 1; ibiny <= covShapeMuKin->GetNbinsY(); ++ibiny)
1048  covMatShapeMuKin[ibinx - 1][ibiny - 1] = covShapeMuKin->GetBinContent(ibinx, ibiny);
1049  for (int ibinx = 1; ibinx <= covShapeENu->GetNbinsX(); ++ibinx)
1050  for (int ibiny = 1; ibiny <= covShapeENu->GetNbinsY(); ++ibiny)
1051  covMatShapeENu[ibinx - 1][ibiny - 1] = covShapeENu->GetBinContent(ibinx, ibiny);
1052  for (int ibinx = 1; ibinx <= covShapeQ2->GetNbinsX(); ++ibinx)
1053  for (int ibiny = 1; ibiny <= covShapeQ2->GetNbinsY(); ++ibiny)
1054  covMatShapeQ2[ibinx - 1][ibiny - 1] = covShapeQ2->GetBinContent(ibinx, ibiny);
1055  TMatrixD covMatShapeMuKinFull = covMatShapeMuKin + statMatMuKin;
1056  TMatrixD covMatShapeENuFull = covMatShapeENu + statMatENu;
1057  TMatrixD covMatShapeQ2Full = covMatShapeQ2 + statMatQ2;
1058  TMatrixD invCovMatShapeMuKin = covMatShapeMuKinFull;
1059  TMatrixD invCovMatShapeENu = covMatShapeENuFull;
1060  TMatrixD invCovMatShapeQ2 = covMatShapeQ2Full;
1061  invCovMatShapeMuKin.Invert();
1062  invCovMatShapeENu .Invert();
1063  invCovMatShapeQ2 .Invert();
1064 
1065  // Plot and save covariances
1066  TDirectory * covarianceDir = outFile->mkdir("covariance");
1067  covarianceDir->cd();
1068  PlotCovMatrices(covMatTotalENuFull, invCovMatTotalENu, outputFolder, "enuTotal", "E_{#nu} Total");
1069  PlotCovMatrices(covMatTotalQ2Full, invCovMatTotalQ2, outputFolder, "q2Total", "q^{2} Total");
1070  PlotCovMatrices(covMatTotalMuKinFull, invCovMatTotalMuKin, outputFolder, "mukinTotal", "Muon Kinematic Total");
1071  PlotCovMatrices(covMatShapeENuFull, invCovMatShapeENu, outputFolder, "enuShape", "E_{#nu} Shape-only");
1072  PlotCovMatrices(covMatShapeQ2Full, invCovMatShapeQ2, outputFolder, "q2Shape", "q^{2} Shape-only");
1073  PlotCovMatrices(covMatShapeMuKinFull, invCovMatShapeMuKin, outputFolder, "mukinShape", "Muon Kinematic Shape-only");
1074 
1075  // Set nominal error bars to the sqrt of covariance diagonals
1076  SetErrorFromCovDiag(nomMuKinShape1D , covMatShapeMuKinFull);
1077  SetErrorFromCovDiag(nomENuShape_nobininfo, covMatShapeENuFull);
1078  SetErrorFromCovDiag(nomQ2Shape_nobininfo , covMatShapeQ2Full);
1079 
1080  // Load the right binnings from the nominal
1081  TH2 * nomMuKinShape = Construct2D(nomMuKinShape1D);
1082  TH1 * nomENuShape = ConstructENu(nomENuShape_nobininfo);
1083  TH1 * nomQ2Shape = ConstructQ2(nomQ2Shape_nobininfo);
1084  SigmaDivideENu(nomENuShape); // ENu then gets scaled.
1085 
1086  // // Load a second copy with the stats errors included
1087  // TH1 * nomMuKinShapeStats1D = (TH1*) nomMuKinShape1D ->Clone("hmukin_cv_2D_shape_stats");
1088  // TH1 * nomENuShapeStats_nobininfo = (TH1*) nomENuShape_nobininfo->Clone("henu_cv_enuaxis_shape_stats");
1089  // TH1 * nomQ2ShapeStats_nobininfo = (TH1*) nomQ2Shape_nobininfo ->Clone("hq2_cv_q2axis_shape_stats");
1090  // SetErrorFromCovDiag(nomMuKinShapeStats1D , covMatShapeMuKinStats);
1091  // SetErrorFromCovDiag(nomENuShapeStats_nobininfo, covMatShapeENuStats);
1092  // SetErrorFromCovDiag(nomQ2ShapeStats_nobininfo , covMatShapeQ2Stats);
1093  // TH2 * nomMuKinShapeStats = Construct2D(nomMuKinShapeStats1D);
1094  // TH1 * nomENuShapeStats = ConstructENu(nomENuShapeStats_nobininfo);
1095  // TH1 * nomQ2ShapeStats = ConstructQ2(nomQ2ShapeStats_nobininfo);
1096  // SigmaDivideENu(nomENuShapeStats); // ENu then gets scaled.
1097 
1098  // Set titles, axis, markers, etc.
1099  nomMuKinShape->SetTitle((";" + muAngTitle + ";" + mukeTitle + ";" + xsec2DTitle).c_str());
1100  nomENuShape-> SetTitle((";" + enuTitle + ";" + xsecENuTitle).c_str());
1101  nomQ2Shape-> SetTitle((";" + q2Title + ";" + xsecQ2Title).c_str());
1102  nomMuKinShape->SetMarkerStyle(20);
1103  nomENuShape-> SetMarkerStyle(20);
1104  nomQ2Shape-> SetMarkerStyle(20);
1105 
1106  // Save to the output
1107  TDirectory * shapeUncertResults = xsecDir->mkdir("shapeUncertResults");
1108  shapeUncertResults->cd();
1109  nomMuKinShape->Write("MuKin");
1110  nomENuShape->Write("ENu");
1111  nomQ2Shape->Write("Q2");
1112 
1113  /////////////////////////////////
1114  // Make map of covariance bin to global bin
1115  /////////////////////////////////
1116  TH1I * covTo2DBin = MakeCovTo2DBin(nomMuKinTotal, covMatTotalMuKin.GetNrows());
1117  outFile->cd();
1118  covTo2DBin->Write("CovTo2DGlobalBin");
1119  bin2DToCov->Write("Bin2DToCov");
1120 
1121  /////////////////////////////////
1122  // Load the Generators
1123  /////////////////////////////////
1124 
1125  cout << "Started loading generators" << endl;
1126 
1127  // Load Generator Predictions
1128  TFile * gibuuFile = TFile::Open("GiBUU_numuCC_predictions.root");
1129  TFile * neutFile = TFile::Open("NEUT54_numuCC_predictions.root");
1130  TFile * nuwroFile = TFile::Open("NuWro_numuCC_predictions.root");
1131  TFile * genieFile = TFile::Open("GeniePrediction.root");
1132  TFile * genie3File = TFile::Open("Genie3Predictions.root");
1133  assert(gibuuFile && gibuuFile->IsOpen());
1134  assert(neutFile && neutFile->IsOpen());
1135  assert(nuwroFile && nuwroFile->IsOpen());
1136  assert(genieFile && genieFile->IsOpen());
1137  assert(genie3File && genie3File->IsOpen());
1138  map<string, map<string, TH1* > > generators;
1139  generators["GiBUU 2019"] = LoadGenerator(gibuuFile, "gibuu", kBlue, 1E39);
1140  generators["NEUT 5.4.0"] = LoadGenerator(neutFile, "neut", kGreen, 10);
1141  generators["NuWro 2019"] = LoadGenerator(nuwroFile, "nuwro", kViolet, 1E39);
1142  generators["GENIE 3.00.06"] = LoadGeniePrediction(genie3File, "genie3", kOrange, "ppfx");
1143  generators["GENIE 2.12.2 - Untuned"] = LoadGeniePrediction(genieFile, "novatune", kCyan, "ppfx");
1144  generators["GENIE 2.12.2 - NOvA Tune"] = LoadGeniePrediction(genieFile, "untuned", kRed, "both");
1145 
1146  // Copy of each generator that has been normalized to area.
1147  map<string, map<string, TH1* > > generatorsAreaNorm;
1148  for(pair<string, map<string, TH1* > > generator : generators)
1149  for(pair<string, TH1*> var : generator.second)
1150  generatorsAreaNorm[generator.first + " (area norm.)"][var.first] = (TH1*) var.second->Clone((var.second->GetName() + string("_norm")).c_str());
1151  for (pair<string, map<string, TH1* > > generator : generatorsAreaNorm)
1152  {
1153  AreaNormalize2D((TH2*) generator.second["MuKin"], nomMuKinShape, generator.first + " 2D");
1154  AreaNormalize(generator.second["ENu"], nomENuShape, enuBinEdges, generator.first + " ENu");
1155  AreaNormalize(generator.second["Q2"], nomQ2Shape, q2BinEdges, generator.first + " Q2");
1156  }
1157 
1158  // Rebinned generators
1159  map<string, map<string, TH1* > > generatorsVarBin;
1160  generatorsVarBin["GiBUU 2019"] = LoadGenerator(gibuuFile, "gibuu", kBlue, 1E39, "rebin");
1161  generatorsVarBin["NEUT 5.4.0"] = LoadGenerator(neutFile, "neut", kGreen, 10, "rebin");
1162  generatorsVarBin["NuWro 2019"] = LoadGenerator(nuwroFile, "nuwro", kViolet, 1E39, "rebin");
1163  generatorsVarBin["GENIE 3.00.06"] = LoadGeniePrediction(genie3File, "genie3_rebin", kOrange, "ppfx");
1164  generatorsVarBin["GENIE 2.12.2 - Untuned"] = LoadGeniePrediction(genieFile, "novatune_rebin", kCyan, "ppfx"); // Genie are identical
1165  generatorsVarBin["GENIE 2.12.2 - NOvA Tune"] = LoadGeniePrediction(genieFile, "untuned_rebin", kRed, "both");
1166  // Need to hardcode fix a few things
1167  generatorsVarBin["GiBUU 2019"]["ENu"]->Rebin(5);
1168  generatorsVarBin["GiBUU 2019"]["ENu"]->Scale(1./5);
1169  generatorsVarBin["GiBUU 2019"]["ENu"]->GetXaxis()->SetRangeUser(0.0, 5.0);
1170  generatorsVarBin["NuWro 2019"]["ENu"]->Rebin(5);
1171  generatorsVarBin["NuWro 2019"]["ENu"]->Scale(1./5);
1172  // generatorsVarBin["NuWro 2019"]["Q2"]->Scale(1.0/22);
1173  // generatorsVarBin["NEUT 5.4.0"]["Q2"]->Scale(1.0/30);
1174  generatorsVarBin["NEUT 5.4.0"]["ENu"]->Rebin(10);
1175  generatorsVarBin["NEUT 5.4.0"]["ENu"]->Scale(1./10);
1176 
1177  // Copy of each generator that has been normalized to area.
1178  map<string, map<string, TH1* > > generatorsAreaNormVarBin;
1179  for(pair<string, map<string, TH1* > > generator : generatorsVarBin)
1180  for(pair<string, TH1*> var : generator.second)
1181  generatorsAreaNormVarBin[generator.first + " (area norm.)"][var.first] = (TH1*) var.second->Clone((var.second->GetName() + string("_norm")).c_str());
1182  for (pair<string, map<string, TH1* > > generator : generatorsAreaNormVarBin)
1183  {
1184  AreaNormalize2D((TH2*) generator.second["MuKin"], nomMuKinShape, generator.first + " 2D Rebin");
1185  AreaNormalize(generator.second["ENu"], nomENuShape, enuBinEdges, generator.first + " ENu Rebin");
1186  AreaNormalize(generator.second["Q2"], nomQ2Shape, q2BinEdges, generator.first + " Q2 Rebin");
1187  }
1188  c1->SetRightMargin(0.15);
1189 
1190  // Rebinned generators
1191  generatorDir->cd();
1192  map<string, map<string, TH1* > > generatorsRebin = RebinGenerators(generatorsVarBin);
1193  map<string, map<string, TH1* > > generatorsAreaNormRebin = RebinGenerators(generatorsAreaNormVarBin);
1194  generatorDir->mkdir("standard")->cd();
1195  for (auto generator : generators)
1196  for (auto var : generator.second)
1197  var.second->Write((generator.first + "_" + var.first).c_str());
1198  generatorDir->mkdir("normalized")->cd();
1199  for (auto generator : generatorsAreaNorm)
1200  for (auto var : generator.second)
1201  var.second->Write((generator.first + "_" + var.first + "_norm").c_str());
1202  generatorDir->mkdir("databinned")->cd();
1203  for (auto generator : generatorsRebin)
1204  for (auto var : generator.second){
1205  var.second->Write((generator.first + "_" + var.first + "_databins").c_str());
1206  }
1207  generatorDir->mkdir("databinned_normalized")->cd();
1208  for (auto generator : generatorsAreaNormRebin)
1209  for (auto var : generator.second){
1210  var.second->Write((generator.first + "_" + var.first + "_databins_norm").c_str());
1211  }
1212  // Delete the ENu from Chi2
1213  // generatorsRebin["NuWro 2019"]->Delete("ENu");
1214 
1215  /////////////////////////////////
1216  // Calculate Chi2s
1217  /////////////////////////////////
1218  map<string, map<string, double > > generatorChi2s;
1219  map<string, map<string, TH1* > > chi2_min_scaled_generators;
1220  map<string, map<string, double > > generatorScaledChi2s;
1221  map<string, string> gen_short_forms
1222  {
1223  {"GiBUU 2019", "GiBUU"},
1224  {"NEUT 5.4.0", "NEUT"},
1225  {"NuWro 2019", "NuWro"},
1226  {"GENIE 2.12.2 - Untuned", "Genie_Untuned"},
1227  {"GENIE 2.12.2 - NOvA Tune", "Genie_Tuned"},
1228  {"GENIE 3.00.06", "Genie_3"}
1229  };
1230  // for (auto generatorPair : generatorsRebin){
1231  // double enuChi2 = CalculateENuChi2(invCovMatTotalENu, nomENuTotal, generatorPair.second["ENu"], gen_short_forms[generatorPair.first]);
1232  // generatorChi2s[generatorPair.first]["ENu"] = enuChi2;
1233  // c1->SaveAs((outputFolder + "/enu_chi2_contribs_" + gen_short_forms[generatorPair.first] + ".pdf").c_str());
1234 
1235  // // TMinuit * enu_minimize = new TMinuit(1);
1236  // // chi2_min_enu_invCovMat = new TMatrixD(invCovMatTotalENu);
1237  // // chi2_min_enu_data = nomENuTotal;
1238  // // chi2_min_enu_gen = generatorPair.second["ENu"];
1239  // // enu_minimize->SetFCN(fcn_enu);
1240  // // enu_minimize->DefineParameter(0, "scale", 1.0, 10.0, 0., 0.);
1241  // // enu_minimize->Migrad();
1242  // // double enu_scale, enu_scale_err;
1243  // // enu_minimize->GetParameter(0, enu_scale, enu_scale_err);
1244  // // cout << gen_short_forms[generatorPair.first] + " ENu scaling: " << enu_scale << " (error = " << enu_scale_err << ")" << endl;
1245  // // TH1 * scaled_enu_gen = (TH1*)generatorPair.second["ENu"]->Clone((generatorPair.second["ENu"]->GetName() + string("_minchi2")).c_str());
1246  // // scaled_enu_gen->Scale(enu_scale);
1247  // // chi2_min_scaled_generators[generatorPair.first]["ENu"] = scaled_enu_gen;
1248  // // generatorScaledChi2s[generatorPair.first]["ENu"] = CalculateENuChi2(invCovMatTotalENu, nomENuTotal, scaled_enu_gen, gen_short_forms[generatorPair.first] + "_chi2min");
1249  // // c1->SaveAs((outputFolder + "/enu_chi2_contribs_" + gen_short_forms[generatorPair.first] + "_minchi2.pdf").c_str());
1250 
1251  // double q2Chi2 = CalculateQ2Chi2(invCovMatTotalQ2, nomQ2Total, generatorPair.second["Q2"], gen_short_forms[generatorPair.first]);
1252  // generatorChi2s[generatorPair.first]["Q2"] = q2Chi2;
1253  // c1->SaveAs((outputFolder + "/q2_chi2_contribs_" + gen_short_forms[generatorPair.first] + ".pdf").c_str());
1254 
1255  // // TMinuit * q2_minimize = new TMinuit(1);
1256  // // chi2_min_q2_invCovMat = new TMatrixD(invCovMatTotalQ2);
1257  // // chi2_min_q2_data = nomQ2Total;
1258  // // chi2_min_q2_gen = generatorPair.second["Q2"];
1259  // // q2_minimize->SetFCN(fcn_q2);
1260  // // q2_minimize->DefineParameter(0, "scale", 1.0, 10.0, 0., 0.);
1261  // // q2_minimize->Migrad();
1262  // // double q2_scale, q2_scale_err;
1263  // // q2_minimize->GetParameter(0, q2_scale, q2_scale_err);
1264  // // cout << gen_short_forms[generatorPair.first] + " Q2 scaling: " << q2_scale << " (error = " << q2_scale_err << ")" << endl;
1265  // // TH1 * scaled_q2_gen = (TH1*)generatorPair.second["Q2"]->Clone((generatorPair.second["Q2"]->GetName() + string("_minchi2")).c_str());
1266  // // scaled_q2_gen->Scale(q2_scale);
1267  // // chi2_min_scaled_generators[generatorPair.first]["Q2"] = scaled_q2_gen;
1268  // // generatorScaledChi2s[generatorPair.first]["Q2"] = CalculateQ2Chi2(invCovMatTotalQ2, nomQ2Total, scaled_q2_gen, gen_short_forms[generatorPair.first] + "_chi2min");
1269  // // c1->SaveAs((outputFolder + "/q2_chi2_contribs_" + gen_short_forms[generatorPair.first] + "_minchi2.pdf").c_str());
1270 
1271  // cout << "2D Chi2: " << generatorPair.first << endl;
1272  // double mukinChi2 = Calculate2DChi2(invCovMatTotalMuKin, nomMuKinTotal, (TH2*) generatorPair.second["MuKin"], covTo2DBin, gen_short_forms[generatorPair.first]);
1273  // cout << "Found 2D chi2 = " << mukinChi2 << " for generator: " << generatorPair.first << endl;
1274  // generatorChi2s[generatorPair.first]["MuKin"] = mukinChi2;
1275  // c1->SaveAs((outputFolder + "/mukin_chi2_contribs_" + gen_short_forms[generatorPair.first] + ".pdf").c_str());
1276  // }
1277  // map<string, map<string, double > > generatorNormChi2s;
1278  // for (auto generatorPair : generatorsAreaNormRebin){
1279  // string gen_short_form = gen_short_forms[generatorPair.first.substr(0, generatorPair.first.size() - 13)] + "_areanorm";
1280 
1281  // cout << "ENu Chi2: " << generatorPair.first << endl;
1282  // double enuChi2 = CalculateENuChi2(invCovMatShapeENu, nomENuShape, generatorPair.second["ENu"], gen_short_form);
1283  // cout << "Found ENu chi2 = " << enuChi2 << " for generator: " << generatorPair.first << endl;
1284  // generatorNormChi2s[generatorPair.first]["ENu"] = enuChi2;
1285  // c1->SaveAs((outputFolder + "/enu_chi2_contribs_" + gen_short_form + ".pdf").c_str());
1286 
1287  // cout << "Q2 Chi2: " << generatorPair.first << endl;
1288  // double q2chi2 = CalculateQ2Chi2(invCovMatShapeQ2, nomQ2Shape, generatorPair.second["Q2"], gen_short_form);
1289  // cout << "Found Q2 chi2 = " << q2chi2 << " for generator: " << generatorPair.first << endl;
1290  // generatorNormChi2s[generatorPair.first]["Q2"] = q2chi2;
1291  // c1->SaveAs((outputFolder + "/q2_chi2_contribs_" + gen_short_form + ".pdf").c_str());
1292 
1293  // cout << "2D Chi2: " << generatorPair.first << endl;
1294  // double mukinChi2 = Calculate2DChi2(invCovMatShapeMuKin, nomMuKinShape, (TH2*) generatorPair.second["MuKin"], covTo2DBin, gen_short_form);
1295  // cout << "Found 2D chi2 = " << mukinChi2 << " for generator: " << generatorPair.first << endl;
1296  // generatorNormChi2s[generatorPair.first]["MuKin"] = mukinChi2;
1297  // c1->SaveAs((outputFolder + "/mukin_chi2_contribs_" + gen_short_form + ".pdf").c_str());
1298  // }
1299 
1300  /////////////////////////////////
1301  // Make plots
1302  /////////////////////////////////
1303 
1304  string dataTotalTitle = "Data - Total error";
1305  string dataShapeTitle = "Data - Shape only errs";
1306 
1307  // Canvas
1308  TLegend * leg;
1309 
1310  // Draw the full 2D XSec to check
1311  c1->SetRightMargin(0.15);
1312  nomMuKinTotal->Draw("COLZ");
1313  c1->SaveAs((outputFolder + "/xsec_2dProj.pdf").c_str());
1314 
1315  // Tweak the margins again
1316  c1->SetLeftMargin(0.13);
1317  c1->SetRightMargin(0.05);
1318  c1->SetBottomMargin(0.10);
1319 
1320  // 2D Errors
1321  PrintAllSlices(nomMuKinTotal, nomMuKinStats, dataTotalTitle, outputFolder + "/muKinSlices", true, generators);
1322  PrintAllSlices(nomMuKinShape, nomMuKinStats, dataShapeTitle, outputFolder + "/muKinSlices_norm", true, generatorsAreaNorm);
1323  PrintAllSlices(nomMuKinTotal, nomMuKinStats, dataTotalTitle, outputFolder + "/muKinSlices_rebin", false, generatorsRebin, true);
1324  PrintAllSlices(nomMuKinShape, nomMuKinStats, dataShapeTitle, outputFolder + "/muKinSlices_norm_rebin", false, generatorsAreaNormRebin, true);
1325 
1326  // ENu Plots
1327  nomENuTotal->GetYaxis()->SetRangeUser(0.0, nomENuTotal->GetMaximum() * 1.35);
1328  PlotDataGenerators(nomENuTotal, nomENuStats, dataTotalTitle, generators, "ENu", true);
1329  c1->SaveAs((outputFolder + "/results_enu.pdf").c_str());
1330  PlotDataGenerators(nomENuShape, nomENuStats, dataShapeTitle, generatorsAreaNorm, "ENu", true);
1331  c1->SaveAs((outputFolder + "/results_enu_norm.pdf").c_str());
1332  PlotDataGenerators(nomENuTotal, nomENuStats, dataTotalTitle, generatorsRebin, "ENu", false);
1333  c1->SaveAs((outputFolder + "/results_enu_rebin.pdf").c_str());
1334  PlotDataGenerators(nomENuShape, nomENuStats, dataShapeTitle, generatorsAreaNormRebin, "ENu", false);
1335  c1->SaveAs((outputFolder + "/results_enu_norm_rebin.pdf").c_str());
1336 
1337  // Q2 Plots
1338  nomQ2Total->GetYaxis()->SetRangeUser(0.0, nomQ2Total->GetMaximum() * 1.25);
1339  nomQ2Shape->GetYaxis()->SetRangeUser(0.0, nomQ2Shape->GetMaximum() * 1.25);
1340  PlotDataGenerators(nomQ2Total, nomQ2Stats, dataTotalTitle, generators, "Q2", true);
1341  c1->SaveAs((outputFolder + "/results_q2.pdf").c_str());
1342  PlotDataGenerators(nomQ2Shape, nomQ2Stats, dataShapeTitle, generatorsAreaNorm, "Q2", true);
1343  c1->SaveAs((outputFolder + "/results_q2_norm.pdf").c_str());
1344  PlotDataGenerators(nomQ2Total, nomQ2Stats, dataTotalTitle, generatorsRebin, "Q2", false);
1345  c1->SaveAs((outputFolder + "/results_q2_rebin.pdf").c_str());
1346  PlotDataGenerators(nomQ2Shape, nomQ2Stats, dataShapeTitle, generatorsAreaNormRebin, "Q2", false);
1347  c1->SaveAs((outputFolder + "/results_q2_norm_rebin.pdf").c_str());
1348 
1349  // Ratios
1350  c1->Clear();
1351  delete c1;
1352  c1 = new TCanvas("c1", "c1", 1600, 1320);
1353  leg = new TLegend();
1354  leg->AddEntry(nomENuTotal, dataTotalTitle.c_str());
1355  PlotGenRatios(generatorsRebin, nomENuTotal, nomENuStats, "ENu", &leg);
1356  c1->SaveAs((outputFolder + "/results_enu_ratios.pdf").c_str());
1357 
1358  c1->Clear();
1359  leg->Clear();
1360  delete leg;
1361  leg = new TLegend();
1362  leg->AddEntry(nomENuShape, dataShapeTitle.c_str());
1363  PlotGenRatios(generatorsAreaNormRebin, nomENuShape, nomENuStats, "ENu", &leg);
1364  c1->SaveAs((outputFolder + "/results_enu_shape_ratios.pdf").c_str());
1365 
1366  c1->Clear();
1367  leg->Clear();
1368  delete leg;
1369  leg = new TLegend();
1370  leg->AddEntry(nomQ2Total, dataTotalTitle.c_str());
1371  PlotGenRatios(generatorsRebin, nomQ2Total, nomQ2Stats, "Q2", &leg);
1372  c1->SaveAs((outputFolder + "/results_q2_ratios.pdf").c_str());
1373 
1374  c1->Clear();
1375  leg->Clear();
1376  delete leg;
1377  leg = new TLegend();
1378  leg->AddEntry(nomQ2Shape, dataShapeTitle.c_str());
1379  PlotGenRatios(generatorsAreaNormRebin, nomQ2Shape, nomQ2Stats, "Q2", &leg);
1380  c1->SaveAs((outputFolder + "/results_q2_shape_ratios.pdf").c_str());
1381 
1382  outFile->Close();
1383 }
void fcn_enu(int &npar, double *gin, double &f, double *par, int iflag)
pair< TH1 *, TH1 * > GeneratorRatio(TH1 *gen, TH1 *data)
enum BeamMode kOrange
TH1 * SetErrorFromCovDiag(TH1 *hist, TH2 *mat)
TMatrixD * chi2_min_enu_invCovMat
enum BeamMode kRed
map< int, int > MuTperCosBins
TH1 * chi2_min_enu_gen
TMatrixD * chi2_min_q2_invCovMat
correl_xv GetYaxis() -> SetDecimals()
const Var weight
string mukeTitle
map< string, TH1 * > LoadGenerator(TFile *file, string label, Color_t col, double scaling=-1, string histLabel="default")
const std::vector< double > angBinEdges_reduced
T sqrt(T number)
Definition: d0nt_math.hpp:156
TH1 * ConstructQ2(TH1 *histIn)
TH1 * ratio(TH1 *h1, TH1 *h2)
void OverlayGenerators(map< string, map< string, TH1 * > > generators, string var, bool drawSpline=false, TLegend **leg=NULL, double binLowEdge=-1.)
constexpr T pow(T x)
Definition: pow.h:75
int Get2DBinArea(TH2 *hist, int binxLow, int binxHigh, int binyLow, int binyHigh)
void PlotDataGenerators(TH1 *dataSystPlusStat, TH1 *dataStat, string dataTitle, map< string, map< string, TH1 * > > &generatorMap, string varName, bool drawSpline=false)
map< string, map< string, string > > generator_hist_names
Int_t par
Definition: SimpleIterate.C:24
TCanvas * canv
OStream cerr
Definition: OStream.cxx:7
void abs(TH1 *hist)
string filename
Definition: shutoffs.py:106
const std::vector< double > mukeBinEdges_reduced
int GetBinByEdge(TAxis *axis, double binEdge, EBinEdge findBinBy=kBinLowEdge)
const std::vector< double > enuBinEdges
map< string, TH1 * > LoadGeniePrediction(TFile *file, string label, Color_t col, string weight)
fVtxDx SetMarkerStyle(20)
gargamelle SetTitle("Gargamelle #nu_{e} CC data")
double chi2()
int GetCovBin(TH1I *covTo2DBin, int globalBin)
void SigmaDivideENu(TH1 *hist)
void scaling(TH1D *hIn, const double shape_scale)
double CalculateENuChi2(const TMatrixD &invCovMat, TH1 *data, TH1 *generator, string label)
const char * label
const XML_Char const XML_Char * data
Definition: expat.h:268
const int nbins
Definition: cellShifts.C:15
string enuTitle
std::string GetName(int i)
void PrintAllSlices(TH2 *histSystPlusStat, TH2 *histStat, string histLabel, string filename, bool drawSpline, map< string, map< string, TH1 * > > generators={}, bool printRatios=false)
TFile * outFile
Definition: PlotXSec.C:135
TH2 * Construct2D(TH1 *histIn)
TH1 * chi2_min_q2_gen
Float_t E
Definition: plot.C:20
TMatrixD StatMatFromHist(TH1 *hist)
TH2I * Make2DToCovBin()
Int_t col[ntarg]
Definition: Style.C:29
correl_xv GetXaxis() -> SetDecimals()
const std::vector< double > q2BinEdges
hmean SetLineWidth(2)
const std::vector< double > mukeBinEdges
void fcn_q2(int &npar, double *gin, double &f, double *par, int iflag)
const std::vector< double > angBinEdges
TH1 * ConstructENu(TH1 *histIn)
void AreaNormalize(TH1 *generator, TH1 *data, const vector< double > &binning, const string label)
const double j
Definition: BetheBloch.cxx:29
const std::string covarianceDir
string xsec2DTitle
string xsecENuTitle
OStream cout
Definition: OStream.cxx:6
void PlotCovariance(string outputFolder=".")
double Calculate1DChi2_simple(const TMatrixD &invCovMat, const TH1 *data, const TH1 *generator)
map< string, map< string, TH1 * > > RebinGenerators(map< string, map< string, TH1 * > > generators)
TH1 * chi2_min_enu_data
enum BeamMode kViolet
EBinEdge
double epsilon
string cv_name
A script which takes covariance matrices, and constructs our final pretty plots for NumuCC...
string q2Title
double Calculate2DChi2(const TMatrixD &invCovMat, TH2 *data, TH2 *generator, TH1 *covToGlobalBinMap, string label)
TH1I * MakeCovTo2DBin(TH2 *hist, int nbins)
const hit & b
Definition: hits.cxx:21
double CalculateQ2Chi2(const TMatrixD &invCovMat, TH1 *data, TH1 *generator, string label)
TFile * file
Definition: cellShifts.C:17
assert(nhit_max >=nhit_nbins)
string muAngTitle
c1
Definition: demo5.py:24
TLegend * AutoPlaceLegend(double dx, double dy, double yPin)
Create a legend, maximizing distance from all histograms.
Definition: Plots.cxx:715
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
void PlotCovMatrices(TMatrixD &cov, TMatrixD &inv, string outputFolder, string label, string title)
TH1 * chi2_min_q2_data
enum BeamMode kGreen
enum BeamMode kBlue
def entry(str)
Definition: HTMLTools.py:26
void PlotGenRatios(map< string, map< string, TH1 * > > generators, TH1 *dataFull, TH1 *dataStat, string var, TLegend **leg=NULL, double binLowEdge=-1.)
void AreaNormalize2D(TH2 *generator, TH2 *data, const string label)
Eigen::MatrixXd mat
string xsecQ2Title
fvar< T > inv(const fvar< T > &x)
Definition: inv.hpp:12
enum BeamMode string