CovMxSurface.cxx
Go to the documentation of this file.
1 #include <iomanip>
2 
3 #include "TDirectory.h"
4 #include "TKey.h"
5 #include "TH2.h"
6 #include "TCanvas.h"
7 #include "TGraph.h"
8 #include "TMultiGraph.h"
9 #include "TLegend.h"
10 #include "TStyle.h"
11 #include "TROOT.h"
12 
13 #include "CAFAna/Core/Progress.h"
14 #include "CAFAna/Core/Utilities.h"
15 #include "CAFAna/Core/IFitVar.h"
16 #include "CAFAna/Core/ISyst.h"
21 
22 #include "OscLib/OscCalcSterile.h"
23 
24 namespace ana {
25 
26  //---------------------------------------------------------------------------
27  /// Default constructor
29  osc::OscCalcSterile* calc, const IFitVar* xvar,
30  const IFitVar* yvar, std::vector<const IFitVar*> fit_vars,
31  std::map<const IFitVar*, std::vector<double> > fit_seeds,
32  std::vector<const ISyst*> prof_systs, SystShifts syst_seeds,
33  bool save_spectra, const ISyst* pull_syst, std::vector<IPrediction*> preds,
34  std::vector<Spectrum*> data, std::vector<Spectrum*> cosmic, size_t patch,
35  size_t n_patches) :
36  fSurface(nullptr), fFitTime(nullptr) {
37 
38  // Create surfaces
39  fSurface = MakeTH2D(UniqueName().c_str(), "Fit surface", xbins, ybins);
40  fFitTime = MakeTH2D(UniqueName().c_str(), "Fit time", xbins, ybins);
41  for (auto var : fit_vars)
42  fFitVarSurface.push_back(MakeTH2D(UniqueName().c_str(), var->ShortName().c_str(),
43  xbins, ybins));
44  for (auto syst : prof_systs)
45  fSystPullSurface.push_back(MakeTH2D(UniqueName().c_str(), syst->ShortName().c_str(),
46  xbins, ybins));
47 
48  // Length of longest fit var name, for debug output formatting
49  size_t pad_width = (*std::max_element(fit_vars.begin(), fit_vars.end(),
50  [] (const IFitVar* lhs, const IFitVar* rhs) {
51  return lhs->ShortName().length() < rhs->ShortName().length();
52  }))->ShortName().length();
53 
55  for (auto var : fit_vars) {
56  if (fit_seeds.count(var) && !fit_seeds[var].empty()) {
57  std::cout << "Seed values for " << std::left << std::setw(pad_width)
58  << var->ShortName() << ": ";
59  for (double seed : fit_seeds[var])
60  std::cout << seed << " ";
62  }
63  }
65 
66  // Figure out which parts of the surface to run over
67  size_t total = xbins.NBins() * ybins.NBins();
68  double stride = double(total)/double(n_patches);
69 
70  // If slice and nslice are not set, we run over everything
71  size_t first, last;
72  if (patch == 999 && n_patches == 999) {
73  first = 0;
74  last = total;
75  }
76 
77  // Otherwise, we figure out the first and last like this
78  else {
79  first = size_t(stride * patch); // this is the first point to run
80  last = size_t(stride * (patch+1)); // this is the count AFTER the last point to run
81  }
82 
83  Progress* p = nullptr;
84  if (!RunningOnGrid()) p = new Progress("Producing surface");
85 
86  for (size_t i = 0; i < total; ++i) {
87 
88  if (i < first) continue;
89  if (i >= last) break;
90  size_t x = (i / ybins.NBins()) + 1;
91  size_t y = (i % ybins.NBins()) + 1;
92  double xval = fSurface->GetXaxis()->GetBinCenter(x);
93  double yval = fSurface->GetYaxis()->GetBinCenter(y);
94 
95  if (RunningOnGrid()) {
97  std::cout << " *********************************************" << std::endl;
98  std::cout << " Bin X Index: " << x
99  << " Bin Y Index: " << x
100  << " Bin Center X: " << xval << " "
101  << " Bin Center Y: " << yval << std::endl;
102  }
103 
104  // Set variables for this bin
105  xvar->SetValue(calc, xval);
106  yvar->SetValue(calc, yval);
107 
108  // Are we saving spectra here?
109  std::vector<std::tuple<std::string, SystShifts, osc::IOscCalcAdjustable*>> calcs;
110  if (save_spectra) calcs.push_back(std::make_tuple("Before fit",
111  SystShifts(syst_seeds), new osc::OscCalcSterile(*calc)));
112 
113  // Do the fit, save the output
114  MinuitFitter fitter(expt, fit_vars, prof_systs);
115  time_t before, after;
116  time(&before);
117  double chi = fitter.Fit(calc, syst_seeds, fit_seeds, {}, IFitter::kQuiet)->EvalMetricVal();
118  time(&after);
119 
120  if (save_spectra) {
121  calcs.push_back(std::make_tuple("After fit", SystShifts(syst_seeds),
122  new osc::OscCalcSterile(*calc)));
123  DrawSurfacePoint(preds, calcs, data, cosmic, x, y);
124  }
125 
126  fSurface->SetBinContent(x, y, chi);
127  double seconds = difftime(after, before);
128  fFitTime->SetBinContent(x, y, seconds);
129  for (size_t it = 0; it < fit_vars.size(); ++it)
130  fFitVarSurface[it]->SetBinContent(x, y, fit_vars[it]->GetValue(calc));
131  for (size_t it = 0; it < prof_systs.size(); ++it)
132  fSystPullSurface[it]->SetBinContent(x, y, syst_seeds.GetShift(prof_systs[it]));
133 
134  for (size_t it = 0; it < prof_systs.size(); ++it)
135  fSystPullSurface[it]->SetBinContent(x, y, syst_seeds.GetShift(prof_systs[it]));
136 
137  if (RunningOnGrid()) {
138  for (auto var : fit_vars)
139  std::cout << " Fit value for " << std::left << std::setw(pad_width)
140  << var->ShortName() << " : " << var->GetValue(calc) << std::endl;
141  std::cout << " ChiSq from fit: " << chi << std::endl;
142  }
143  else p->SetProgress(double(i - first) / double(last - first));
144 
145  } // for surface point i
146 
147  if (!RunningOnGrid()) p->Done();
148 
149  } // CovMxSurface constructor
150 
151  //---------------------------------------------------------------------------
152  /// Default destructor
154 
155  delete fSurface;
156  delete fFitTime;
157 
158  for (TH2D* var : fFitVarSurface) delete var;
159 
160  for (TH2D* syst : fSystPullSurface) delete syst;
161 
162  for (auto it : fSurfacePoints)
163  if (it.second) delete it.second;
164 
165  } // CovMxSurface destructor
166 
167  //---------------------------------------------------------------------------
168  /// Get contours from surface
169  std::vector<std::vector<TGraph*> > CovMxSurface::GetContours(
170  std::vector<double> up) {
171 
172  std::vector<std::vector<TGraph*> > ret(up.size(), std::vector<TGraph*>());
173 
174  TCanvas c(UniqueName().c_str(), "");
175 
176  // Loop over each up-value and get the contours
177 
178  for (size_t i = 0; i < up.size(); ++i) {
179 
180  std::vector<TGraph*> g;
181 
182  // This is the best way to get contours from a surface, because
183  // ROOT is a nightmare
184 
185  TH2D* h = (TH2D*)fSurface->Clone();
186  double min = fSurface->GetBinContent(fSurface->GetMinimumBin());
187  h->SetContour(1);
188  h->SetContourLevel(0, min + TMath::ChisquareQuantile(up[i], 1));
189  h->Draw("cont list same");
190  c.Update();
191  TObjArray* contour_array =
192  (TObjArray*)gROOT->GetListOfSpecials()->FindObject("contours");
193  TList* contour_list = (TList*)contour_array->At(0);
194  TIter next(contour_list);
195  TGraph* graph = nullptr;
196  while ( (graph = (TGraph*)next()) )
197  g.push_back((TGraph*)graph->Clone());
198  ret[i] = g;
199  delete h;
200 
201  } // for up value
202 
203  return ret;
204 
205  } // function CovMxSurface::GetContours
206 
207  //---------------------------------------------------------------------------
208  /// Draw contours from multiple surfaces on the same canvas
209  TCanvas* CovMxSurface::DrawContours(std::vector<CovMxSurface*> surfs,
210  std::vector<double> up, std::vector<std::string> names) {
211 
212  TCanvas* ret = new TCanvas(UniqueName().c_str(), "", 1200, 675);
213 
214  // Define colours for different surfaces
215 
216  std::vector<int> colours = { kAzure+2, kRed+2, kGreen+2, kBlack };
217  if (surfs.size() < 1 || surfs.size() > 4)
218  assert(false && "You must provide 1-4 surfaces.");
219 
220  // Define line styles for different up-values
221 
222  std::vector<int> styles = { kSolid, kDotted, kDashed };
223  if (up.size() < 1 || up.size() > 3)
224  assert(false && "You must provide 1-3 up-values.");
225 
226  // Check the labels make sense if provided
227 
228  if (!names.empty() && names.size() != surfs.size())
229  assert(false && "List of labels does not match list of surfaces!");
230 
231  // Draw a template histogram
232 
233  TH2D* h = (TH2D*)surfs[0]->fSurface->Clone();
234  h->Draw("axis");
235 
236  // Draw legend
237 
238  TLegend* l = new TLegend(0.2, 0.65, 0.4, 0.85);
239 
240  // Draw graphs
241 
242  for (size_t it_surf = 0; it_surf < surfs.size(); ++it_surf) {
243  std::vector<std::vector<TGraph*> > graphs
244  = surfs[it_surf]->GetContours(up);
245  for (size_t it_up = 0; it_up < up.size(); ++it_up) {
246 
247  for (TGraph* g : graphs[it_up]) {
248  g->SetLineWidth(3);
249  g->SetLineColor(colours[it_surf]);
250  g->SetLineStyle(styles[it_up]);
251  g->Draw("c same");
252  }
253  if (!graphs[it_up].empty() && !names.empty())
254  l->AddEntry(graphs[it_up][0], Form("%s (%.0f%% CL)",
255  names[it_surf].c_str(), 100*up[it_up]), "l");
256 
257  } // for up-value
258  } // for surface
259 
260  l->SetFillStyle(0);
261  l->SetTextFont(42);
262  l->SetBorderSize(0);
263  l->SetTextSize(0.03);
264  l->Draw("same");
265 
266  ret->SetLogx(true);
267  ret->SetLogy(true);
268 
269  return ret;
270 
271  } // function CovMxSurface::DrawContours
272 
273  //---------------------------------------------------------------------------
274  /// Standard SaveTo implementation
275  void CovMxSurface::SaveTo(TDirectory* dir) const {
276 
277  TObjString tag("CovMxSurface");
278  dir->WriteTObject(&tag, "type");
279  dir->WriteTObject(fSurface, "Surface");
280  dir->WriteTObject(fFitTime, "Fit time");
281 
282  TDirectory* var_dir = dir->mkdir("Fit vars");
283  for (auto surf : fFitVarSurface)
284  var_dir->WriteTObject(surf, surf->GetTitle());
285 
286  TDirectory* syst_dir = dir->mkdir("Syst pulls");
287  for (auto surf : fSystPullSurface)
288  syst_dir->WriteTObject(surf, surf->GetTitle());
289 
290  TDirectory* spec_dir = dir->mkdir("Spectra");
291  for (auto spectrum : fSurfacePoints)
292  spec_dir->WriteTObject(spectrum.second, spectrum.first.c_str());
293 
294  } // function CovMxSurface::SaveTo
295 
296  //---------------------------------------------------------------------------
297  /// Standard LoadFrom implementation
298  std::unique_ptr<CovMxSurface> CovMxSurface::LoadFrom(TDirectory* dir) {
299 
300  DontAddDirectory guard;
301 
302  TObjString* tag = (TObjString*)dir->Get("type");
303  assert(tag and tag->GetString() == "CovMxSurface");
304 
305  TH2D* surface = (TH2D*)dir->Get("Surface");
306  assert(surface);
307 
308  TH2D* time = (TH2D*)dir->Get("Fit time");
309  assert(time);
310 
311  std::vector<TH2D*> vars;
312  TDirectory* var_dir = dir->GetDirectory("Fit vars");
313  if (var_dir) {
314  TIter it(var_dir->GetListOfKeys());
315  TKey* key = nullptr;
316  while ((key = (TKey*)it())) {
317  TH2D* h_tmp = (TH2D*)key->ReadObj();
318  assert(h_tmp);
319  vars.push_back(h_tmp);
320  }
321  }
322 
323  std::vector<TH2D*> pulls;
324  TDirectory* pull_dir = dir->GetDirectory("Syst pulls");
325  if (pull_dir) {
326  TIter it(pull_dir->GetListOfKeys());
327  TKey* key = nullptr;
328  while ((key = (TKey*)it())) {
329  TH2D* h_tmp = (TH2D*)key->ReadObj();
330  assert(h_tmp);
331  pulls.push_back(h_tmp);
332  }
333  }
334 
335  std::map<std::string, TCanvas*> specs;
336  TDirectory* spec_dir = dir->GetDirectory("Spectra");
337  if (spec_dir) {
338  TIter it(spec_dir->GetListOfKeys());
339  TKey* key = nullptr;
340  while ((key = (TKey*)it())) {
341  std::string name = key->GetName();
342  TCanvas* c = (TCanvas*)key->ReadObj();
343  assert(c);
344  specs[name] = c;
345  }
346  }
347 
348  std::unique_ptr<CovMxSurface> ret(new CovMxSurface(surface, vars,
349  pulls, time, specs));
350 
351  return std::move(ret);
352 
353  } // function CovMxSurface::LoadFrom
354 
355  //---------------------------------------------------------------------------
356  /// Private constructor used by LoadFrom
357  CovMxSurface::CovMxSurface(TH2D* surface, std::vector<TH2D*> fit_vars,
358  std::vector<TH2D*> syst_pulls, TH2D* fit_time,
359  std::map<std::string, TCanvas*> spectra) :
360  fSurface(surface), fFitVarSurface(fit_vars), fSystPullSurface(syst_pulls),
361  fFitTime(fit_time), fSurfacePoints(spectra)
362  {}
363 
364  //---------------------------------------------------------------------------
365  /// Draw spectra pre- and post-fit
366  void CovMxSurface::DrawSurfacePoint(std::vector<IPrediction*> preds,
367  std::vector<std::tuple<std::string, SystShifts,
368  osc::IOscCalcAdjustable*>> calcs, std::vector<Spectrum*> data,
369  std::vector<Spectrum*> cosmic, size_t x, size_t y) {
370 
371  DontAddDirectory guard;
372 
373  std::vector<int> colours = { kAzure-2, kRed+2, kGreen+3, kMagenta+2 };
374  int cosmic_colour = kGray+2;
375 
376  // Get the number of detectors and universes
377 
378  if (calcs.empty()) throw std::runtime_error("Spectrum vector is empty.");
379  size_t n_universes = calcs.size();
380  if (n_universes > colours.size())
381  throw std::runtime_error("Only " + std::to_string(colours.size())
382  + " universes per plot are supported!");
383 
384  size_t n_detectors = preds.size();
385 
386  // Use OscCalcs and simulation to get spectra
387 
388  std::vector<std::vector<Spectrum>> specs(n_universes);
389  for (size_t i = 0; i < n_universes; ++i) {
390  for (auto pred : preds) {
391  specs[i].push_back(pred->PredictSyst(std::get<2>(calcs[i]),
392  std::get<1>(calcs[i])));
393  }
394  }
395 
396  // Create data and cosmic histograms
397  std::vector<TH1D*> h_data;
398  std::vector<TH1D*> h_cosmic;
399  for (size_t i = 0; i < data.size(); ++i) {
400  h_data.push_back(data[i]->ToTH1(data[i]->POT()));
401  if (cosmic[i]) h_cosmic.push_back(cosmic[i]->ToTH1(data[i]->Livetime(), kLivetime));
402  else h_cosmic.push_back(nullptr);
403  }
404 
405  // Define shape of canvas
406  TCanvas* c = new TCanvas(UniqueName().c_str(), "c", 500 * n_detectors, 800);
407  c->Divide(n_detectors, 3);
408 
409  // Pad margins
410  std::vector<float> margin = { 0.15, 0.05 };
411  float ts = 20, ls = 14;
412 
413  for (size_t d = 0; d < n_detectors; ++d)
414  {
415  // Get pad x ranges for this detector
416  float xmin = float(d) / float(n_detectors);
417  float xmax = float(d+1) / float(n_detectors);
418 
419  // Set size of spectrum pad
420  c->cd(d+1);
421  gPad->SetPad(xmin, 0.45, xmax, 1);
422  gPad->SetLeftMargin(margin[0]);
423  gPad->SetRightMargin(margin[1]);
424 
425  // Set size of ratio pad
426  c->cd(n_detectors+d+1);
427  gPad->SetPad(xmin, 0.25, xmax, 0.45);
428  gPad->SetLeftMargin(margin[0]);
429  gPad->SetRightMargin(margin[1]);
430 
431  // Set size of oscillation pad
432  c->cd((2*n_detectors)+d+1);
433  gPad->SetPad(xmin, 0.05, xmax, 0.25);
434  gPad->SetLeftMargin(margin[0]);
435  gPad->SetRightMargin(margin[1]);
436 
437  } // for detector i
438 
439  TLegend* l = new TLegend(0.5, 0.65, 0.85, 0.85);
440  l->SetFillStyle(0);
441 
442  // Create histogram vectors and set their sizes
443 
444  std::vector<std::vector<TH1D*>> hists, ratio_hists;
445 
446  hists.resize(n_universes);
447  for (std::vector<TH1D*>& h : hists)
448  h.resize(n_detectors);
449  ratio_hists.resize(n_universes);
450  for (std::vector<TH1D*>& h : ratio_hists)
451  h.resize(n_detectors);
452  std::vector<TMultiGraph*> graphs(n_detectors, nullptr);
453 
454  // Y axis max for each detector's plot
455  std::vector<double> ymax, rmax;
456  std::vector<bool> rescaled;
457 
458  // Make data ratio histograms
459 
460  std::vector<TH1D*> data_ratio;
461  data_ratio.resize(n_detectors);
462  for (size_t d = 0; d < n_detectors; ++d)
463  {
464  data_ratio[d] = (TH1D*)h_data[d]->Clone();
465  for (int b = 1; b <= data_ratio[d]->GetNbinsX(); ++b)
466  {
467  double val = h_data[d]->GetBinContent(b);
468  double err = h_data[d]->GetBinError(b);
469  data_ratio[d]->SetBinContent(b, 1);
470 
471  if (val == 0)
472  data_ratio[d]->SetBinError(b, 0);
473  else
474  data_ratio[d]->SetBinError(b, err / val);
475 
476  } // for bin b
477  } // for detector d
478 
479  // Loop over each detector
480 
481  std::vector<bool> is_nd(n_detectors);
482 
483  for (size_t d = 0; d < n_detectors; ++d)
484  {
485  // Loop over each universe, get histograms and get ymax
486 
487  std::vector<double> m, r;
488 
489  // Rescale data histogram
490  is_nd[d] = (h_data[d]->GetMaximum() > 1e4);
491  if (is_nd[d]) h_data[d]->Scale(1e-4);
492 
493  for (size_t u = 0; u < n_universes; ++u)
494  {
495  hists[u][d] = (specs[u][d].ToTH1(data[d]->POT()));
496 
497  // If it's a near detector spectrum, rescale it
498  if (is_nd[d]) hists[u][d]->Scale(1e-4);
499  m.push_back(hists[u][d]->GetMaximum());
500 
501  // Clone to create ratio histogram
502  ratio_hists[u][d] = (TH1D*)hists[u][d]->Clone();
503  TH1D* h = ratio_hists[u][d];
504  for (int b = 1; b <= h->GetNbinsX(); ++b)
505  {
506  double val;
507  if (h_cosmic[d])
508  val = (hists[u][d]->GetBinContent(b) + h_cosmic[d]->GetBinContent(b)) / h_data[d]->GetBinContent(b);
509  else
510  val = hists[u][d]->GetBinContent(b) / h_data[d]->GetBinContent(b);
511  if (std::isinf(val) || std::isnan(val)) {
512  h->SetBinContent(b, 1);
513  h->SetBinError(b, 0);
514  }
515  else {
516  h->SetBinContent(b, val);
517  r.push_back(fabs(1-val));
518  }
519 
520  } // for bin b
521 
522  } // for universe u
523 
524  // Get maximum data point
525 
526  int maxdata = h_data[d]->GetMaximumBin();
527  m.push_back(h_data[d]->GetBinContent(maxdata) + h_data[d]->GetBinError(maxdata));
528 
529  // Account for data errors in ratio plot range
530 
531  // for (int b = 1; b <= data_ratio[d]->GetNbinsX(); ++b)
532  // r.push_back(fabs(data_ratio[d]->GetBinError(b)));
533 
534  ymax.push_back(*std::max_element(m.begin(), m.end()));
535  rmax.push_back(*std::max_element(r.begin(), r.end()));
536 
537  // Initialise oscillation calculator multigraphs
538 
539  graphs[d] = new TMultiGraph();
540 
541  } // for detector d
542 
543  // Loop over each detector to draw the cosmics
544  for (size_t d = 0; d < n_detectors; ++d) {
545 
546  // Draw the cosmic spectrum
547  c->cd(d+1);
548  if (h_cosmic[d]) {
549 
550  TH1D* hc = h_cosmic[d];
551  hc->SetLineColor(cosmic_colour);
552  hc->SetLineWidth(2);
553  hc->SetMinimum(0);
554  hc->SetMaximum(1.1 * ymax[d]);
555  hc->Draw("hist same");
556 
557  // Also need to set the axis labels and titles here!
558  if (is_nd[d])
559  hc->SetTitle(";Energy deposited in detector (GeV);10^{4} events");
560  else
561  hc->SetTitle(";Energy deposited in detector (GeV);Events");
562 
563  hc->GetXaxis()->SetTitleFont(43);
564  hc->GetYaxis()->SetTitleFont(43);
565  hc->GetXaxis()->SetLabelFont(43);
566  hc->GetYaxis()->SetLabelFont(43);
567 
568  hc->GetXaxis()->SetTitleSize(ts);
569  hc->GetYaxis()->SetTitleSize(ts);
570  hc->GetXaxis()->SetLabelSize(ls);
571  hc->GetYaxis()->SetLabelSize(ls);
572 
573  hc->GetXaxis()->SetTitleOffset(1.5);
574  hc->GetYaxis()->SetTitleOffset(2.5);
575  }
576  } // for detector d
577 
578  // Loop over each universe
579  for (size_t u = 0; u < n_universes; ++u) {
580  // Loop over each detector
581  for (size_t d = 0; d < n_detectors; ++d) {
582  // Format and draw the spectrum
583 
584  c->cd(d+1);
585  TH1D* h = hists[u][d];
586  if (h_cosmic[d]) h->Add(h_cosmic[d]);
587  h->SetLineColor(colours[u]);
588  h->SetLineWidth(2);
589  h->SetMinimum(0);
590  h->SetMaximum(1.1 * ymax[d]);
591 
592  h->Draw("hist same");
593 
594  // Now draw the ratio plot
595  h = ratio_hists[u][d];
596  c->cd(n_detectors+d+1);
597  h->SetLineColor(colours[u]);
598  h->SetLineWidth(2);
599  h->SetMinimum(1 - (1.1 * rmax[d]));
600  h->SetMaximum(1 + (1.1 * rmax[d]));
601 
602  h->SetTitle(";;Ratio");
603  h->GetYaxis()->SetTitleOffset(2.5);
604 
605  h->Draw("hist same");
606 
607  h->GetXaxis()->SetTitleFont(43);
608  h->GetYaxis()->SetTitleFont(43);
609  h->GetXaxis()->SetLabelFont(43);
610  h->GetYaxis()->SetLabelFont(43);
611 
612  h->GetXaxis()->SetTitleSize(ts);
613  h->GetYaxis()->SetTitleSize(ts);
614  h->GetXaxis()->SetLabelSize(ls);
615  h->GetYaxis()->SetLabelSize(ls);
616 
617  c->Update();
618 
619  // Now the oscillation plot
620  float xmin = h_data[d]->GetXaxis()->GetXmin();
621  float xmax = h_data[d]->GetXaxis()->GetXmax();
622  auto calc = std::get<2>(calcs[u]);
623  std::vector<float> x, y;
624  double l = calc->GetL();
625  if (is_nd[d]) calc->SetL(1);
626  else calc->SetL(810);
627  for (size_t k = 0; k < 1000; ++k)
628  {
629  x.push_back(xmin + (0.001 * (xmax - xmin) * (float(k) + 0.5)));
630  y.push_back(calc->P(14, 0, x.back()));
631  }
632  calc->SetL(l);
633 
634  TGraph* g = new TGraph(x.size(), &x[0], &y[0]);
635  g->SetLineWidth(2);
636  g->SetLineColor(colours[u]);
637  graphs[d]->Add(g);
638 
639  } // for detector d
640 
641  l->AddEntry(hists[u][0], std::get<0>(calcs[u]).c_str(), "l");
642 
643  } // for universe u
644 
645  // Cosmic legend entry
646  for (auto h : h_cosmic) {
647  if (h) {
648  l->AddEntry(h, "Cosmic background", "l");
649  break;
650  }
651  }
652 
653  if (h_data.size() != n_detectors)
654  throw std::runtime_error("Number of data hists does not match number of samples!");
655 
656  // Loop and draw the data and oscillation graph for each sample
657  for (size_t d = 0; d < n_detectors; ++d)
658  {
659  c->cd(d+1);
660  h_data[d]->SetMarkerSize(2);
661  h_data[d]->Draw("e0 same");
662 
663  c->cd((2*n_detectors)+d+1);
664  graphs[d]->Draw("al");
665 
666  graphs[d]->SetTitle(";;1 - P(#nu_{#mu} #rightarrow #nu_{s})");
667  graphs[d]->GetYaxis()->SetTitleOffset(2.5);
668  graphs[d]->GetXaxis()->SetRangeUser(h_data[d]->GetXaxis()->GetXmin(),
669  h_data[d]->GetXaxis()->GetXmax());
670 
671  graphs[d]->GetXaxis()->SetTitleFont(43);
672  graphs[d]->GetYaxis()->SetTitleFont(43);
673  graphs[d]->GetXaxis()->SetLabelFont(43);
674  graphs[d]->GetYaxis()->SetLabelFont(43);
675 
676  graphs[d]->GetXaxis()->SetTitleSize(ts);
677  graphs[d]->GetYaxis()->SetTitleSize(ts);
678  graphs[d]->GetXaxis()->SetLabelSize(ls);
679  graphs[d]->GetYaxis()->SetLabelSize(ls);
680 
681  } // for detector d
682 
683  l->AddEntry(h_data[0], "Data", "lep");
684  c->cd(1);
685  l->Draw();
686 
687  c->Update();
688 
689  std::ostringstream name;
690  name << "spectrum_" << x << "_" << y;
691  fSurfacePoints[name.str()] = c;
692 
693  } // function DrawSurfacePoint
694 
695  //---------------------------------------------------------------------------
696  /// Take a concatenated fake data spectrum and split it up
697 /* std::vector<TH1D*> CovMxSurface::SplitFakeData(TH1D* h,
698  std::vector<Binning> b)
699  {
700  // Define split-up histograms
701 
702  std::vector<TH1D*> h_split;
703  size_t n_bins = 0;
704  for (size_t i = 0; i < b.size(); ++i) {
705  h_split.push_back(HistCache::New(UniqueName(), b[i]));
706  n_bins += b[i].NBins();
707  }
708 
709  // Check binning is consistent
710 
711  if (h->GetNbinsX() != int(n_bins))
712  throw std::runtime_error("There are " + std::to_string(h->GetNbinsX())
713  + " bins in merged histogram but a total of " + std::to_string(n_bins)
714  + " bins across the binning vector provided.");
715 
716  int bin_count = 0;
717 
718  for (size_t i = 0; i < b.size(); ++i) {
719  for (int j = 1; j <= h_split[i]->GetNbinsX(); ++j) {
720  double val = h->GetBinContent(bin_count+j);
721  h_split[i]->SetBinContent(j, val);
722  h_split[i]->SetBinError(j, sqrt(val));
723  }
724  bin_count += h_split[i]->GetNbinsX();
725  }
726 
727  return h_split;
728 
729  } // function SplitFakeData
730 */
731 } // namespace ana
~CovMxSurface()
Default destructor.
CovMxSurface(Binning xbins, Binning ybins, IExperiment *expt, osc::OscCalcSterile *calc, const IFitVar *xvar, const IFitVar *yvar, std::vector< const IFitVar * > fit_vars, std::map< const IFitVar *, std::vector< double > > fit_seeds, std::vector< const ISyst * > prof_systs, SystShifts syst_seeds, bool save_spectra=false, const ISyst *pull_syst=nullptr, std::vector< IPrediction * > preds={}, std::vector< Spectrum * > data={}, std::vector< Spectrum * > cosmic={}, size_t patch=999, size_t n_patches=999)
Default constructor.
const XML_Char * name
Definition: expat.h:151
int isinf(const stan::math::var &a)
Definition: std_isinf.hpp:16
def after(args)
Definition: fabricate.py:1416
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::map< std::string, double > xmax
set< int >::iterator it
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
void DrawSurfacePoint(std::vector< IPrediction * > preds, std::vector< std::tuple< std::string, SystShifts, osc::IOscCalcAdjustable * > > calcs, std::vector< Spectrum * > data, std::vector< Spectrum * > cosmic, size_t x, size_t y)
Draw spectra pre- and post-fit.
bool RunningOnGrid()
Is this a grid (condor) job?
Definition: Utilities.cxx:363
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
const char * p
Definition: xmltok.h:285
Adapt the PMNS_Sterile calculator to standard interface.
std::vector< TH2D * > fSystPullSurface
Fit variable best fit points.
Definition: CovMxSurface.h:68
TString hists[nhists]
Definition: bdt_com.C:3
virtual double P(int flavBefore, int flavAfter, double E) override
E in GeV; flavors as PDG codes (so, neg==>antinu)
Definition: OscCalcDumb.h:21
osc::OscCalcDumb calc
::xsd::cxx::tree::time< char, simple_type > time
Definition: Database.h:194
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
void SaveTo(TDirectory *dir) const
Standard SaveTo implementation.
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
const XML_Char const XML_Char * data
Definition: expat.h:268
std::vector< std::vector< TGraph * > > GetContours(std::vector< double > up)
Get contours from surface.
Double_t ymax
Definition: plot.C:25
expt
Definition: demo5.py:34
Sum up livetimes from individual cosmic triggers.
const int xbins
Definition: MakePlots.C:82
unsigned int seed
Definition: runWimpSim.h:102
correl_xv GetXaxis() -> SetDecimals()
virtual std::unique_ptr< IFitSummary > Fit(osc::IOscCalcAdjustable *seed, SystShifts &bestSysts=junkShifts, const SeedList &seedPts=SeedList(), const std::vector< SystShifts > &systSeedPts={}, Verbosity verb=kVerbose) const
Master fitting method. Depends on FitHelper and FitHelperSeeded.
Definition: IFitter.cxx:69
T GetShift(const ISyst *syst) const
Float_t d
Definition: plot.C:236
const std::map< std::pair< std::string, std::string >, Variable > vars
OStream cout
Definition: OStream.cxx:6
std::vector< double > POT
cout<< "--"<< endl;for(Int_t iP=1;iP<=hyz->GetNbinsX();iP++){for(Int_t iC=1;iC<=hyz->GetNbinsY();iC++){if(hyv->GetBinContent(iP, iC)>-999){goal_hyv-> SetBinContent(iP, iC,-(dy[iP-1][iC-1]))
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
const int ybins
TH2D * MakeTH2D(const char *name, const char *title, const Binning &binsx, const Binning &binsy)
Utility function to avoid 4-way combinatorial explosion on the bin types.
Definition: UtilsExt.cxx:86
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:104
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void SetProgress(double frac)
Update the progress fraction between zero and one.
Definition: Progress.cxx:41
const std::string & ShortName() const
Definition: IFitVar.h:36
TDirectory * dir
Definition: macro.C:5
virtual void SetValue(osc::IOscCalcAdjustable *osc, double val) const =0
std::vector< TH2D * > fFitVarSurface
Surface histogram.
Definition: CovMxSurface.h:67
int NBins() const
Definition: Binning.h:29
Base class defining interface for experiments.
Definition: IExperiment.h:14
const hit & b
Definition: hits.cxx:21
std::map< std::string, TCanvas * > fSurfacePoints
Surface tracking fit time for each surface point.
Definition: CovMxSurface.h:72
assert(nhit_max >=nhit_nbins)
A simple ascii-art progress bar.
Definition: Progress.h:9
TRandom3 r(0)
Interface definition for fittable variables.
Definition: IFitVar.h:16
def ls(target="")
Definition: g4zmq.py:69
def save_spectra(fname, spectra, groups)
Definition: core.py:116
static std::unique_ptr< CovMxSurface > LoadFrom(TDirectory *dir)
Standard LoadFrom implementation.
TH2D * fFitTime
Systematic pulls.
Definition: CovMxSurface.h:70
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
Float_t e
Definition: plot.C:35
void Done()
Call this when action is completed.
Definition: Progress.cxx:90
void next()
Definition: show_event.C:84
Definition: styles.py:1
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
static TCanvas * DrawContours(std::vector< CovMxSurface * > surfs, std::vector< double > up, std::vector< std::string > names={})
Draw contours from multiple surfaces on the same canvas.
float h_data[xbins][ybins]
surf
Definition: demo5.py:35
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17