Plots.cxx
Go to the documentation of this file.
2 
4 
5 #include "CAFAna/Core/Ratio.h"
6 #include "CAFAna/Core/Spectrum.h"
8 #include "Utilities/func/MathUtil.h"
9 
10 #include "TCanvas.h"
11 #include "TFeldmanCousins.h"
12 #include "TGraph.h"
13 #include "TGraphAsymmErrors.h"
14 #include "TH1.h"
15 #include "TH2.h"
16 #include "THStack.h"
17 #include "THLimitsFinder.h"
18 #include "TLegend.h"
19 #include "TLine.h"
20 #include "TList.h"
21 #include "TMarker.h"
22 #include "TPaveText.h"
23 #include "TText.h"
24 
25 #include <algorithm>
26 #include <iostream>
27 #include <cmath>
28 #include <ctime>
29 #include <functional>
30 #include <memory>
31 
32 namespace ana
33 {
34  //----------------------------------------------------------------------
35  TH1* DataMCComparison(const Spectrum& data, const Spectrum& mc, EBinType bintype)
36  {
37  TH1* ret = 0;
38 
39  TH1* hMC = mc.ToTH1(data.POT(), kTotalMCColor, kSolid, kPOT, bintype);
40 
41  TH1* hData = data.ToTH1(data.POT(), kBlack, kSolid, kPOT, bintype);
42  hData->Sumw2();
43  hData->SetMarkerStyle(kFullCircle);
44 
45  // Need to draw the biggest thing first to get correct axis limits
46  if(hMC->GetMaximum() > hData->GetMaximum()+sqrt(hData->GetMaximum())){
47  ret = hMC;
48  hMC->Draw("hist");
49  hData->Draw("ep same");
50  }
51  else{
52  ret = hData;
53  hData->Draw("ep");
54  hMC->Draw("hist same");
55  hData->Draw("ep same"); // Data always goes on top
56  }
57 
58  gPad->Update();
59 
60  return ret;
61  }
62 
63  //----------------------------------------------------------------------
65  const IPrediction* mc,
67  const SystShifts & shifts,
68  EBinType bintype)
69  {
70  return DataMCComparison(data, shifts.IsNominal() ? mc->Predict(calc) : mc->PredictSyst(calc, shifts), bintype);
71  }
72 
73  //----------------------------------------------------------------------
75  {
76  TH1* ret = 0;
77 
78  TH1* hMC = mc.ToTH1(mc.POT());
79  hMC->SetLineColor(kTotalMCColor);
80 
81  TH1* hData = data.ToTH1(data.POT());
82  hData->Sumw2();
83  hData->SetMarkerStyle(kFullCircle);
84 
85  hMC->Scale(hData->Integral()/hMC->Integral());
86 
87  // Need to draw the biggest thing first to get correct axis limits
88  if(hMC->GetMaximum() > hData->GetMaximum()+sqrt(hData->GetMaximum())){
89  ret = hMC;
90  hMC->Draw("hist");
91  hData->Draw("ep same");
92  }
93  else{
94  ret = hData;
95  hData->Draw("ep");
96  hMC->Draw("hist same");
97  hData->Draw("ep same"); // Data always goes on top
98  }
99 
100  gPad->Update();
101 
102  return ret;
103  }
104 
105  //----------------------------------------------------------------------
107  const IPrediction* mc,
109  {
110  return DataMCComparisonAreaNormalized(data, mc->Predict(calc));
111  }
112 
113  //----------------------------------------------------------------------
115  const IPrediction* mc,
117  {
118  TH1* ret = 0;
119 
120  TH1* hMC = mc->Predict(calc).ToTH1(data.POT());
121  hMC->SetLineColor(kTotalMCColor);
122 
123  TH1* hNC = mc->PredictComponent(calc,
125  Current::kNC,
126  Sign::kBoth).ToTH1(data.POT());
127  hNC->SetLineColor(kNCBackgroundColor);
128 
129  TH1* hCC = mc->PredictComponent(calc,
131  Current::kCC,
132  Sign::kBoth).ToTH1(data.POT());
133  hCC->SetLineColor(kNumuBackgroundColor);
134 
135  TH1* hBeam = mc->PredictComponent(calc,
137  Current::kCC,
138  Sign::kBoth).ToTH1(data.POT());
139  hBeam->SetLineColor(kBeamNueBackgroundColor);
140 
141  TH1* hData = data.ToTH1(data.POT());
142  hData->SetMarkerStyle(kFullCircle);
143 
144  // Need to draw the biggest thing first to get correct axis limits
145  if(hMC->GetMaximum() > hData->GetMaximum()+sqrt(hData->GetMaximum())){
146  ret = hMC;
147  hMC->Draw("hist");
148  hData->Draw("ep same");
149  }
150  else{
151  ret = hData;
152  hData->Draw("ep");
153  hMC->Draw("hist same");
154  hData->Draw("ep same"); // Data always goes on top
155  }
156 
157  hNC->Draw("hist same");
158  hCC->Draw("hist same");
159  hBeam->Draw("hist same");
160 
161  gPad->Update();
162 
163  return ret;
164  }
165 
166  //----------------------------------------------------------------------
167  void DataMCRatio(const Spectrum& data,
168  const IPrediction* mc,
170  double miny, double maxy)
171  {
172  DataMCRatio(data, mc->Predict(calc), miny, maxy);
173  }
174 
175  //----------------------------------------------------------------------
176  void DataMCRatio(const Spectrum& data,
177  const Spectrum& mc,
178  double miny, double maxy)
179  {
180  Ratio ratio(data, mc);
181  TH1* h = ratio.ToTH1();
182  h->GetYaxis()->SetTitle("Data / MC");
183  h->GetYaxis()->SetRangeUser(miny, maxy);
184  h->Draw();
185 
186  TLine* one = new TLine(h->GetXaxis()->GetXmin(), 1,
187  h->GetXaxis()->GetXmax(), 1);
188  one->SetLineWidth(2);
189  one->SetLineStyle(7);
190  one->Draw();
191 
192  gPad->Update();
193  }
194 
195  //----------------------------------------------------------------------
197  const IPrediction* mc,
199  double miny, double maxy)
200  {
201  DataMCAreaNormalizedRatio(data, mc->Predict(calc), miny, maxy);
202  }
203 
204  //----------------------------------------------------------------------
206  const Spectrum& mc,
207  double miny, double maxy)
208  {
209  Spectrum mcScaled = mc;
210  mcScaled.OverridePOT(mcScaled.POT() * mc.Integral(1) / data.Integral(1));
211 
212  DataMCRatio(data, mcScaled, miny, maxy);
213  }
214 
215  //----------------------------------------------------------------------
216  void RatioPlot(const Spectrum& data,
217  const Spectrum& expected,
218  const Spectrum& fit,
219  double miny, double maxy)
220  {
221  Ratio fitRatio(fit, expected);
222  Ratio dataRatio(data, expected);
223 
224  TH1* h = fitRatio.ToTH1();
225  h->GetYaxis()->SetTitle("Ratio to expectation");
226  h->GetYaxis()->SetRangeUser(miny, maxy);
227  h->SetLineColor(kTotalMCColor);
228  h->Draw("][");
229 
230  h = dataRatio.ToTH1();
231  h->SetMarkerStyle(kFullCircle);
232  h->Draw("ep same");
233 
234  TLine* one = new TLine(h->GetXaxis()->GetXmin(), 1,
235  h->GetXaxis()->GetXmax(), 1);
236  one->SetLineWidth(2);
237  one->SetLineStyle(7);
238  one->Draw();
239 
240  gPad->Update();
241  }
242 
243  //----------------------------------------------------------------------
245  const IPrediction* pred,
246  std::vector<const ISyst*> systs,
247  std::vector<TH1*> &hUps,
248  std::vector<TH1*> &hDns,
249  double pot,
250  Flavors::Flavors_t flav,
253  {
254  for (const ISyst* syst: systs) {
255  SystShifts shifts;
256 
257  shifts.SetShift(syst, +1);
258  TH1* hUp = pred->PredictComponentSyst(calc, shifts, flav, curr, sign).ToTH1(pot);
259 
260  shifts.SetShift(syst, -1);
261  TH1* hDn = pred->PredictComponentSyst(calc, shifts, flav, curr, sign).ToTH1(pot);
262 
263  hUps.push_back(hUp);
264  hDns.push_back(hDn);
265  }
266 
267  return;
268  }
269 
270  //----------------------------------------------------------------------
272  const IPrediction* pred,
273  std::vector<const ISyst*> systs,
274  const SystShifts &bfshifts,
275  std::vector<TH1*> &hUps,
276  std::vector<TH1*> &hDns,
277  double pot,
278  Flavors::Flavors_t flav,
281  {
282  for (const ISyst* syst: systs) {
283  SystShifts shifts;
284 
285  shifts.SetShift(syst, +1);
286  TH1* hUp = pred->PredictComponentSyst(calc, shifts, flav, curr, sign).ToTH1(pot);
287  hUp->Add( pred->PredictComponentSyst(calc, SystShifts::Nominal(), flav, curr, sign).ToTH1(pot), -1 ); // remove the nominal pred
288  hUp->Add( pred->PredictComponentSyst(calc, bfshifts, flav, curr, sign).ToTH1(pot), +1 ); // add on the BF shifted pred
289 
290  shifts.SetShift(syst, -1);
291  TH1* hDn = pred->PredictComponentSyst(calc, shifts, flav, curr, sign).ToTH1(pot);
292  hDn->Add( pred->PredictComponentSyst(calc, SystShifts::Nominal(), flav, curr, sign).ToTH1(pot), -1 ); // remove the nominal pred
293  hDn->Add( pred->PredictComponentSyst(calc, bfshifts, flav, curr, sign).ToTH1(pot), +1 ); // add on the BF shifted pred
294 
295  hUps.push_back(hUp);
296  hDns.push_back(hDn);
297  }
298 
299  return;
300  }
301 
302 
303  //----------------------------------------------------------------------
305  const std::vector<const ISyst*>& systs,
307  double pot,
308  int col, int errCol, float headroom,
309  bool newaxis,
310  EBinType bintype)
311  {
312 
313  Spectrum nom = pred->Predict(calc);
314 
315  std::vector<Spectrum> ups, dns;
316 
317  for(const ISyst* syst: systs){
318  SystShifts shifts;
319  shifts.SetShift(syst, +1);
320  ups.push_back(pred->PredictSyst(calc, shifts));
321  shifts.SetShift(syst, -1);
322  dns.push_back(pred->PredictSyst(calc, shifts));
323  }
324 
325  return PlotWithSystErrorBand(nom, ups, dns, pot, col, errCol, headroom, newaxis, bintype);
326 
327 
328  }
329 
330  //----------------------------------------------------------------------
331 
332  TGraphAsymmErrors* PlotWithSystErrorBand(const Spectrum& nominal,
333  const std::vector<Spectrum>& upShifts,
334  const std::vector<Spectrum>& downShifts,
335  double pot, int col, int errCol,
336  float headroom, bool newaxis,
337  EBinType bintype)
338  {
339 
340  TH1* nom = nominal.ToTH1(pot, kBlack, kSolid, kPOT, bintype);
341 
342  std::vector<TH1*> ups, dns;
343 
344  for(const auto& upShift:upShifts) ups.push_back(upShift.ToTH1(pot, kBlack, kSolid, kPOT, bintype));
345  for(const auto& downShift:downShifts) dns.push_back(downShift.ToTH1(pot, kBlack, kSolid, kPOT, bintype));
346 
347  TGraphAsymmErrors* g = PlotWithSystErrorBand(nom, ups, dns, col, errCol, headroom, newaxis);
348 
349  for(TH1* up: ups) delete up;
350  for(TH1* dn: dns) delete dn;
351 
352  return g;
353  }
354 
355  //----------------------------------------------------------------------
356 
357  TGraphAsymmErrors* PlotWithSystErrorBand(TH1*& nom,
358  std::vector<TH1*>& ups,
359  std::vector<TH1*>& dns,
360  int col, int errCol,
361  float headroom, bool newaxis)
362  {
363  if(col == -1){
364  col = kTotalMCColor;
365  errCol = kTotalMCErrorBandColor;
366  }
367  else if(errCol == -1) errCol = col-7; // hopefully a lighter version
368 
369  nom->SetLineColor(col);
370  nom->GetXaxis()->CenterTitle();
371  nom->GetYaxis()->CenterTitle();
372  if(newaxis) nom->Draw("hist ]["); // Set the axes up
373 
374  double yMax = nom->GetBinContent(nom->GetMaximumBin());
375 
376  TGraphAsymmErrors* g = new TGraphAsymmErrors;
377 
378  for(int binIdx = 0; binIdx < nom->GetNbinsX()+2; ++binIdx){
379  const double y = nom->GetBinContent(binIdx);
380  g->SetPoint(binIdx, nom->GetXaxis()->GetBinCenter(binIdx), y);
381 
382  const double w = nom->GetXaxis()->GetBinWidth(binIdx);
383 
384  double errUp = 0;
385  double errDn = 0;
386 
387  for(unsigned int systIdx = 0; systIdx < ups.size(); ++systIdx){
388  double hi = ups[systIdx]->GetBinContent(binIdx)-y;
389  double lo = dns[systIdx]->GetBinContent(binIdx)-y;
390 
391  // If both shifts are in the same direction use the larger
392  double min = std::min(hi,lo);
393  double max = std::max(hi,lo);
394  if(max < 0) max=0;
395  if(min > 0) min=0;
396 
397  errUp += max*max;
398  errDn += min*min;
399  } // end for systIdx
400 
401  g->SetPointError(binIdx, w/2, w/2, sqrt(errDn), sqrt(errUp));
402  } // end for i
403 
404  g->SetFillColor(errCol);
405  g->Draw("e2 same");
406  g->GetYaxis()->SetRangeUser(0.00001, headroom*yMax);
407  nom->GetYaxis()->SetRangeUser(0.0001, headroom*yMax);
408 
409  nom->Draw("hist same");
410 
411  gPad->RedrawAxis();
412 
413  return g;
414  }
415 
416  //----------------------------------------------------------------------
418  const std::vector<const ISyst*>& systs,
420  double pot,
421  int col, int errCol, float headroom,
422  bool newaxis, EBinType bintype)
423  {
424 
425  Spectrum nom = pred->Predict(calc);
426 
427  std::vector<Spectrum> ups, dns;
428 
429  for(const ISyst* syst: systs){
430  SystShifts shifts;
431  shifts.SetShift(syst, +1);
432  ups.push_back(pred->PredictSyst(calc, shifts));
433  shifts.SetShift(syst, -1);
434  dns.push_back(pred->PredictSyst(calc, shifts));
435  }
436 
437  PlotWithAreaSystErrorBand(nom, ups, dns, pot, col, errCol, headroom, newaxis, bintype);
438 
439 
440  }
441 
442  //----------------------------------------------------------------------
443  void PlotWithAreaSystErrorBand(const Spectrum& nominal,
444  std::vector<Spectrum> upShifts,
445  std::vector<Spectrum> downShifts,
446  double pot, int col, int errCol,
447  float headroom, bool newaxis,
448  EBinType bintype)
449  {
450  double nomIntegral = nominal.Integral(pot);
451  if (nomIntegral > 0)
452  {
453  // ugh so much boilerplate
454  for (auto & collection : std::vector<std::reference_wrapper<std::vector<Spectrum>>>{ upShifts, downShifts })
455  {
456  for (Spectrum & shiftSpec : collection.get())
457  {
458  // this is a dirty trick, but these Spectra are going to be thrown out once this method exits anyway
459  shiftSpec.OverridePOT(shiftSpec.POT() * shiftSpec.Integral(pot) / nomIntegral);
460  }
461  }
462  }
463 
464  PlotWithSystErrorBand(nominal, upShifts, downShifts, pot, col, errCol, headroom, newaxis, bintype);
465  }
466 
467  //----------------------------------------------------------------------
469  const std::vector<Spectrum>& upShifts,
470  const std::vector<Spectrum>& downShifts,
471  const Spectrum& alternative,
472  const std::vector<Spectrum>& upShiftsAlt,
473  const std::vector<Spectrum>& downShiftsAlt,
474  double pot, int col1, int col2,
475  float headroom, bool newaxis,
476  EBinType bintype)
477  {
478  TH1* nom = nominal.ToTH1(pot, kBlack, kSolid, kPOT, bintype);
479  TH1* alt = alternative.ToTH1(pot, kBlack, kSolid, kPOT, bintype);
480  std::vector<TH1*> ups, dns;
481  std::vector<TH1*> upsA, dnsA;
482  for(const auto& upShift:upShifts) ups.push_back(upShift.ToTH1(pot, kBlack, kSolid, kPOT, bintype));
483  for(const auto& downShift:downShifts) dns.push_back(downShift.ToTH1(pot, kBlack, kSolid, kPOT, bintype));
484  for(const auto& upShiftAlt:upShiftsAlt) upsA.push_back(upShiftAlt.ToTH1(pot, kBlack, kSolid, kPOT, bintype));
485  for(const auto& downShiftAlt:downShiftsAlt) dnsA.push_back(downShiftAlt.ToTH1(pot, kBlack, kSolid, kPOT, bintype));
486  nom->SetLineColor(col1);
487  nom->GetXaxis()->CenterTitle();
488  nom->GetYaxis()->CenterTitle();
489  if(newaxis) nom->Draw("hist"); // Set the axes up
490  alt->SetLineColor(col2);
491  double yMax = nom->GetBinContent(nom->GetMaximumBin());
492  TGraphAsymmErrors* g = new TGraphAsymmErrors;
493  TGraphAsymmErrors* g2 = new TGraphAsymmErrors;
494  for(int binIdx = 0; binIdx < nom->GetNbinsX()+2; ++binIdx){
495  const double y = nom->GetBinContent(binIdx);
496  const double y2 = alt->GetBinContent(binIdx);
497  g->SetPoint(binIdx, nom->GetXaxis()->GetBinCenter(binIdx), y);
498  g2->SetPoint(binIdx, nom->GetXaxis()->GetBinCenter(binIdx), y2);
499  const double w = nom->GetXaxis()->GetBinWidth(binIdx);
500  double errUp = 0, errUp2 = 0;
501  double errDn = 0, errDn2 = 0;
502  for(unsigned int systIdx = 0; systIdx < ups.size(); ++systIdx){
503  double hi = ups[systIdx]->GetBinContent(binIdx)-y;
504  double hi2 = upsA[systIdx]->GetBinContent(binIdx)-y2;
505  double lo = y-dns[systIdx]->GetBinContent(binIdx);
506  double lo2 = y2-dnsA[systIdx]->GetBinContent(binIdx);
507  if(lo <= 0 && hi <= 0) std::swap(lo, hi);
508  if(lo2 <= 0 && hi2 <= 0) std::swap(lo2, hi2);
509  errUp += hi*hi;
510  errDn += lo*lo;
511  errUp2 += hi2*hi2;
512  errDn2 += lo2*lo2;
513  // TODO: what happens if they're both high or both low?
514  } // end for systIdx
515  g->SetPointError(binIdx, w/2, w/2, sqrt(errDn), sqrt(errUp));
516  g2->SetPointError(binIdx, w/2, w/2, sqrt(errDn2), sqrt(errUp2));
517  } // end for i
518  g2->SetFillColorAlpha(col2, 0.3);
519  g2->Draw("2");
520  // g2->GetYaxis()->SetRangeUser(0, headroom*yMax);
521  g2->GetYaxis()->SetRangeUser(0.0, 19.0);
522  g->SetFillColorAlpha(col1, 0.3);
523  g->Draw("2");
524  // g->GetYaxis()->SetRangeUser(0, headroom*yMax);
525  // nom->GetYaxis()->SetRangeUser(0, headroom*yMax);
526  g->GetYaxis()->SetRangeUser(0.0, yMax*headroom);
527  nom->GetYaxis()->SetRangeUser(0.0, yMax*headroom);
528  alt->Draw("hist same");
529  nom->Draw("hist same");
530  for(TH1* up: ups) delete up;
531  for(TH1* dn: dns) delete dn;
532  for(TH1* up2: upsA) delete up2;
533  for(TH1* dn2: dnsA) delete dn2;
534  }
535 
536 
537 
538 
539 
540  //----------------------------------------------------------------------
541  TGraphAsymmErrors* PlotWithSystErrorBand_Quant(const int quant,
542  IPrediction* pred,
543  const std::vector<const ISyst*>& systs,
545  double pot,
546  int col, int errCol,
547  float maxy,
548  bool newaxis)
549  {
550 
551  Spectrum nom = pred->Predict(calc);
552 
553  std::vector<Spectrum> ups, dns;
554 
555  for(const ISyst* syst: systs){
556  SystShifts shifts;
557  shifts.SetShift(syst, +1);
558  ups.push_back(pred->PredictSyst(calc, shifts));
559  shifts.SetShift(syst, -1);
560  dns.push_back(pred->PredictSyst(calc, shifts));
561  }
562 
563  return PlotWithSystErrorBand_Quant(quant, nom, ups, dns, pot, col, errCol, maxy, newaxis);
564 
565 
566  }
567 
568  //----------------------------------------------------------------------
569  TGraphAsymmErrors* PlotWithSystErrorBand_Quant(const int quant,
570  const Spectrum& nominal,
571  const std::vector<Spectrum>& upShifts,
572  const std::vector<Spectrum>& downShifts,
573  double pot, int col, int errCol,
574  float maxy, bool newaxis)
575  {
576 
577  TH1* nom = nominal.ToTH1(pot);
578  nom->GetYaxis()->SetTitle("Events/0.1 GeV");
579  std::vector<TH1*> ups, dns;
580 
581  for(const auto& upShift:upShifts) ups.push_back(upShift.ToTH1(pot));
582  for(const auto& downShift:downShifts) dns.push_back(downShift.ToTH1(pot));
583 
584  TGraphAsymmErrors* g = PlotWithSystErrorBand_Quant(quant, nom, ups, dns, col, errCol, maxy, newaxis);
585 
586  for(TH1* up: ups) delete up;
587  for(TH1* dn: dns) delete dn;
588 
589  return g;
590  }
591 
592  //----------------------------------------------------------------------
593  TGraphAsymmErrors* PlotWithSystErrorBand_Quant(const int quant,
594  TH1*& nom,
595  std::vector<TH1*>& ups,
596  std::vector<TH1*>& dns,
597  int col, int errCol,
598  float maxy, bool newaxis)
599  {
600  if(col == -1){
601  col = kTotalMCColor;
602  errCol = kTotalMCErrorBandColor;
603  }
604  else if(errCol == -1) errCol = col-7; // hopefully a lighter version
605 
606  nom->GetYaxis()->SetTitle("Events/0.1 GeV");
607  nom->SetLineColor(col);
608  nom->GetXaxis()->CenterTitle();
609  nom->GetYaxis()->CenterTitle();
610  if(newaxis) nom->Draw("hist ]["); // Set the axes up
611 
612  TGraphAsymmErrors* g = new TGraphAsymmErrors;
613 
614  for(int binIdx = 0; binIdx < nom->GetNbinsX()+2; ++binIdx){
615  const double y = nom->GetBinContent(binIdx);
616  g->SetPoint(binIdx, nom->GetXaxis()->GetBinCenter(binIdx), y);
617 
618  const double w = nom->GetXaxis()->GetBinWidth(binIdx);
619 
620  double errUp = 0;
621  double errDn = 0;
622 
623  for(unsigned int systIdx = 0; systIdx < ups.size(); ++systIdx){
624  double hi = ups[systIdx]->GetBinContent(binIdx)-y;
625  double lo = dns[systIdx]->GetBinContent(binIdx)-y;
626 
627  // If both shifts are in the same direction use the larger
628  double min = std::min(hi,lo);
629  double max = std::max(hi,lo);
630  if(max < 0) max=0;
631  if(min > 0) min=0;
632 
633  errUp += max*max;
634  errDn += min*min;
635  } // end for systIdx
636 
637  g->SetPointError(binIdx, w/2, w/2, sqrt(errDn), sqrt(errUp));
638  } // end for i
639 
640 
641  // setting 2x2 = 4 pads canvas
642  if(quant==1 || quant==2){
643  nom->GetXaxis()->SetLabelSize(0.);
644  }
645  if(quant==2 || quant==4){
646  nom->GetYaxis()->SetLabelSize(0.);
647  }
648  float minx = 0.99 * nom->GetXaxis()->GetXmin();
649  float maxx = 0.99 * nom->GetXaxis()->GetXmax();
650  nom->GetXaxis()->SetLabelOffset(0.01);
651  nom->GetYaxis()->SetLabelOffset(0.01);
652  nom->GetYaxis()->SetTitleOffset(0.50);
653  nom->GetXaxis()->SetTitleSize(0.0);
654  nom->GetYaxis()->SetTitleSize(0.0);
655  nom->GetXaxis()->SetTitleSize(0.0);
656  nom->GetYaxis()->SetTitleSize(0.0);
657  nom->GetYaxis()->SetRangeUser(0.00001, maxy);
658  nom->GetXaxis()->SetRangeUser(minx, maxx);
659  g->GetYaxis() ->SetRangeUser(0.00001, maxy);
660  g->GetXaxis() ->SetRangeUser(minx, maxx);
661  g->SetFillColor(errCol);
662  g->Draw("e2 same");
663  nom->Draw("hist ][ same");
664 
665  return g;
666  }
667 
668 
669 
670 
671 
672 
673  //----------------------------------------------------------------------
674  THStack* ToTHStack(const std::vector<std::pair<const Spectrum&, Color_t>>& s,
675  double pot)
676  {
677  THStack* ret = new THStack;
678  for(auto it: s){
679  TH1* h = it.first.ToTH1(pot);
680  h->SetFillColor(it.second);
681  ret->Add(h, "hist");
682  }
683  return ret;
684  }
685 
686  //----------------------------------------------------------------------
687  /// Helper for \ref AutoPlaceLegend
688  double PointDistanceToBox(double x, double y,
689  double x0, double y0, double x1, double y1)
690  {
691  // Inside
692  if(x > x0 && x < x1 && y > y0 && y < y1) return 0;
693 
694  // Corners
695  double d = util::sqr(x-x0)+util::sqr(y-y0);
696  d = std::min(d, util::sqr(x-x1)+util::sqr(y-y0));
697  d = std::min(d, util::sqr(x-x1)+util::sqr(y-y1));
698  d = std::min(d, util::sqr(x-x0)+util::sqr(y-y1));
699 
700  // Top and bottom edges
701  if(x > x0 && x < x1){
702  d = std::min(d, util::sqr(y-y0));
703  d = std::min(d, util::sqr(y-y1));
704  }
705  // Left and right
706  if(y > y0 && y < y1){
707  d = std::min(d, util::sqr(x-x0));
708  d = std::min(d, util::sqr(x-x1));
709  }
710 
711  return d;
712  }
713 
714  //----------------------------------------------------------------------
715  TLegend* AutoPlaceLegend(double dx, double dy, double yPin)
716  {
717  gPad->Update();
718 
719  // Convert requested width and height into physics coordinates
720  dx *= (gPad->GetX2()-gPad->GetX1());
721  dy *= (gPad->GetY2()-gPad->GetY1());
722 
723  // Range of axes in physics units
724  const double x0 = gPad->GetUxmin();
725  const double x1 = gPad->GetUxmax();
726  const double y0 = gPad->GetUymin();
727  const double y1 = gPad->GetUymax();
728 
729  const double X = x1-x0;
730  const double Y = y1-y0;
731 
732  double bestd = 0;
733  double bestx = 0;
734  double besty = 0;
735  // If we want to pin Y pos, set it to that now.
736  if(yPin >= 0) besty = yPin * (gPad->GetY2()-gPad->GetY1());
737 
738  for(bool fallback: {false, true}){
739  for(double x = x0+dx/2; x < x1-dx/2; x += X/50){
740  for(double y = y0+dy/2; y < y1-dy/2; y += Y/50){
741  double d = 999999;
742 
743  // Repel from edges
744  d = std::min(d, util::sqr((x-dx/2-x0)/X));
745  d = std::min(d, util::sqr((x+dx/2-x1)/X));
746  d = std::min(d, util::sqr((y-dy/2-y0)/Y));
747  d = std::min(d, util::sqr((y+dy/2-y1)/Y));
748  if(d < bestd) continue;
749 
750  TIter next(gPad->GetListOfPrimitives());
751  while(TObject* obj = next()){
752 
753  if(!obj->InheritsFrom(TH1::Class())) continue;
754  if( obj->InheritsFrom(TH2::Class())) continue;
755 
756  TH1* h = (TH1*)obj;
757 
758  for(int n = 1; n <= h->GetNbinsX(); ++n){
759  const double px = h->GetBinCenter(n);
760  const double py = h->GetBinContent(n);
761 
762  if(fallback){
763  d = std::min(d, util::sqr((px-x)/X)+util::sqr((py-y)/Y));
764  }
765  else{
766  d = std::min(d, PointDistanceToBox(px/X, py/Y,
767  (x-dx/2)/X, (y-dy/2)/Y,
768  (x+dx/2)/X, (y+dy/2)/Y));
769  }
770  if(d < bestd) break;
771  }
772  if(d < bestd) break;
773  } // end while
774 
775  if(d > bestd){
776  bestd = d;
777  bestx = x;
778  // Update Y if we're not pinning it.
779  if (yPin < 0) besty = y;
780  }
781  } // end for y
782  } // end for x
783 
784  if(bestd != 0) break; // If we always collide, have to do fallback
785  } // end for fallback
786 
787  // Convert to pad coordinates
788  const double nx = (bestx-gPad->GetX1())/(gPad->GetX2()-gPad->GetX1());
789  const double ny = (besty-gPad->GetY1())/(gPad->GetY2()-gPad->GetY1());
790 
791  const double ndx = dx/(gPad->GetX2()-gPad->GetX1());
792  const double ndy = dy/(gPad->GetY2()-gPad->GetY1());
793 
794  return new TLegend(nx-ndx/2, ny-ndy/2, nx+ndx/2, ny+ndy/2);
795  }
796 
797  //----------------------------------------------------------------------
798  void ErrorBarChart(const std::map<std::string, std::pair<double, double>>& systErrors,
799  const std::pair<double, double> & statErr,
800  const std::string & label )
801  {
802  bool doStat = statErr.first != 0 && statErr.second != 0;
803  int nbins = systErrors.size() + 1; // extra for the "total systematic" bar
804 
805  // Systematics are on top, total systematic and (optional) statistical
806  // error bar are at the bottom.
807  int firstSystBin = 1;
808  if(doStat){
809  ++nbins;
810  ++firstSystBin;
811  }
812 
813  // Sum all the systematics in qaudrature
814  std::pair<double, double> systTot = std::make_pair(0., 0.);
815  for(auto it_systPairs: systErrors)
816  {
817  systTot.first += util::sqr(it_systPairs.second.first);
818  systTot.second += util::sqr(it_systPairs.second.second);
819  }
820  systTot.first = sqrt(systTot.first);
821  systTot.second = sqrt(systTot.second);
822 
823  // x-axis range is a multiple of 5% with at least 2% padding
824  double xrange = 0;
825  while(xrange < std::max(std::max(systTot.first, systTot.second), std::max(statErr.first, statErr.second))+2) xrange += 5;
826 
827  TH2* axes = new TH2F("", (";" + label).c_str(),
828  10, -xrange, +xrange, nbins, 0, nbins);
829  axes->GetXaxis()->CenterTitle();
830 
831  // Put the systematics in a vector so we can sort them by size
832  std::vector<std::pair<std::pair<double, double>, std::string>> systList;
833  for(auto it: systErrors) systList.push_back(std::make_pair(it.second, it.first));
834  std::sort(systList.begin(),
835  systList.end(),
836  [](const auto & errNamePair1, const auto & errNamePair2)
837  {
838  return ( std::max(errNamePair1.first.first, errNamePair1.first.second)
839  < std::max(errNamePair2.first.first, errNamePair2.first.second) );
840  }
841  );
842 
843  // Label the y axis
844  TAxis* yax = axes->GetYaxis();
845  for(unsigned int i = 0; i < systList.size(); ++i)
846  yax->SetBinLabel(i+firstSystBin+1, systList[i].second.c_str());
847 
848  yax->SetBinLabel(firstSystBin, "Total syst. error");
849  if(doStat) yax->SetBinLabel(1, "Statistical error");
850 
851  // Make the labels approximately match the x axis title
852  yax->SetLabelSize(.06);
853 
854  axes->Draw();
855 
856  // Blue boxes for each systematic
857  for(unsigned int i = 0; i < systList.size(); ++i){
858  TBox* box = new TBox(-systList[i].first.first, i+firstSystBin+.2,
859  +systList[i].first.second, i+firstSystBin+.8);
860  std::cout<<"ERROR = (" << systList[i].first.first << ", " << systList[i].first.second << ")" << std::endl;
861  box->SetFillColor(kAzure-8);
862  box->Draw();
863  }
864 
865  // Total systematic error
866  TBox* syst = new TBox(-systTot.first, firstSystBin-.8,
867  +systTot.second, firstSystBin-.2);
868  std::cout<<"TOTAL SYSTEMATIC ERROR = (" << systTot.first << ", " << systTot.second << ")" << std::endl;
869  syst->SetFillColor(kAzure-8);
870  syst->Draw();
871 
872  // Red box for statistical error
873  if(doStat){
874  TBox* stat = new TBox(-statErr.first, .2, +statErr.second, .8);
875  stat->SetFillColor(kRed-7);
876  stat->Draw();
877  }
878 
879  // Separate individual systematics from the total and statistical errors
880  TLine* div = new TLine(-xrange, firstSystBin, +xrange, firstSystBin);
881  div->SetLineWidth(2);
882  div->Draw();
883 
884  // Vertical line marking the symmetry point
885  TLine* zero = new TLine(0, 0, 0, nbins);
886  zero->SetLineStyle(7);
887  zero->SetLineWidth(2);
888  zero->Draw();
889 
890  // Leave enough space for the systematic labels
891  gPad->SetLeftMargin(.25);
892 
893  }
894 
895  //----------------------------------------------------------------------
896  void CountingExperimentErrorBarChart(const std::map<std::string, double>& systbars, double statErr, bool bkgdOrSig, bool shortchart)
897  {
898  std::map<std::string, double> systs = systbars; // breaks shortchart option (shortchart) ? SumNueSecondAnaSysts(systbars) : systbars;
899 
900  std::map<std::string, std::pair<double, double>> systVals;
901  for (const auto & systPair : systbars)
902  systVals[systPair.first] = std::make_pair(systPair.second, systPair.second);
903  auto statVals = std::make_pair(statErr, statErr);
904 
905  std::string label = std::string(bkgdOrSig ? "Signal" : "Background") + " uncertainty (%)";
906  ErrorBarChart(systVals, statVals, label);
907  }
908 
909  //----------------------------------------------------------------------
910  TGraphAsymmErrors* GraphWithPoissonErrors(const TH1* h, bool noErrorsXaxis, bool drawEmptyBins)
911  {
912  TGraphAsymmErrors* gr = new TGraphAsymmErrors(h);
913 
914  TFeldmanCousins fc(0.6827);//1 sigma
915 
916  for(int i = 0; i < h->GetNbinsX(); ++i){
917  double x, y;
918  gr->GetPoint(i, x, y);
919 
920  if ( drawEmptyBins || y!=0 ){
921  if ( y < 50 ) gr->SetPointEYlow(i, y-fc.CalculateLowerLimit(y,0));
922  else gr->SetPointEYlow(i,sqrt(y));
923  if ( y < 30 ) gr->SetPointEYhigh(i,fc.CalculateUpperLimit(y,0)-y);
924  else gr->SetPointEYhigh(i,sqrt(y));
925  }
926 
927  if(noErrorsXaxis){
928  gr->SetPointEXlow(i, 0);
929  gr->SetPointEXhigh(i, 0);
930  } // Do not use bin width as X-error for points
931  }
932  gr->SetMarkerStyle(kFullCircle);
933  gr->SetMarkerColor(1);
934  gr->SetMarkerSize(1);
935  gr->SetLineWidth(2);
936 
937 
938  return gr;
939  }
940 
941  //----------------------------------------------------------------------
942  TGraphAsymmErrors* GraphWithPoissonErrors2(const TH1* h, const TH1* h2, bool noErrorsXaxis, bool drawEmptyBins)
943  {
944  TGraphAsymmErrors* gr = new TGraphAsymmErrors(h);
945  TGraphAsymmErrors* gr2 = new TGraphAsymmErrors();
946 
947  TFeldmanCousins fc(0.6827);//1 sigma
948 
949  int j = 0;
950  for(int i = 0; i < h->GetNbinsX(); ++i){
951  double x, y;
952  gr->GetPoint(i, x, y);
953 
954  y = round(y);
955 
956  if ( (drawEmptyBins && h2->GetBinContent(i+1) > 0.05) || y!=0 ){
957  gr2->SetPoint(j, x, y);
958  gr2->SetPointEXlow(j, gr->GetErrorXlow(i));
959  gr2->SetPointEXhigh(j, gr->GetErrorXhigh(i));
960  if ( y < 50 ) gr2->SetPointEYlow(j, y-fc.CalculateLowerLimit(y,0));
961  else gr2->SetPointEYlow(j,sqrt(y));
962  if ( y < 30 ) gr2->SetPointEYhigh(j,fc.CalculateUpperLimit(y,0)-y);
963  else gr2->SetPointEYhigh(j,sqrt(y));
964  j++;
965  }
966 
967  if(noErrorsXaxis){
968  gr2->SetPointEXlow(i, 0);
969  gr2->SetPointEXhigh(i, 0);
970  } // Do not use bin width as X-error for points
971  }
972 
973  return gr2;
974  }
975 
976 
977  //----------------------------------------------------------------------
978  TGraph* graphAsymmErrorScaled(TH1* histScaled, TH1* hist, double overallScale){
979 
980  for(int bin = 1; bin <= hist->GetNbinsX(); bin++) {
981  if(hist->GetBinContent(bin)==0){
982  hist->SetBinContent(bin,0.001);
983  histScaled->SetBinContent(bin,0.001);
984  }
985  }
986 
987  TFeldmanCousins fc(0.6827); // to get correct poisson error bars
988  TGraphAsymmErrors *datapoisson = new TGraphAsymmErrors(histScaled);
989 
990  for(int bin = 1; bin <= hist->GetNbinsX(); bin++) {
991  //default numu: scale by bin 0.1/width
992  double binWidth = hist->GetBinWidth(bin);
993  double scaleFactor = overallScale / binWidth;
994  double y = hist->GetBinContent(bin);
995 
996  // since this is data it should be a whole number
997  // fc Upper/Lower limit is finicky eg if y = 3.01
998  // the returned value is the same as when y = 4
999  y = round(y);
1000 
1001  double errup = 0, errdn = 0;
1002  if(y < 30) errup = (fc.CalculateUpperLimit(y,0) - y);
1003  else errup = sqrt(y);
1004  if(y < 50) errdn = (y - fc.CalculateLowerLimit(y,0));
1005  else errdn = sqrt(y);
1006 
1007  datapoisson->SetPointEYhigh(bin-1, scaleFactor * errup);
1008  datapoisson->SetPointEYlow (bin-1, scaleFactor * errdn);
1009  }
1010  datapoisson->SetMarkerStyle(kFullCircle);
1011  datapoisson->SetMarkerColor(1);
1012  datapoisson->SetMarkerSize(1);
1013  datapoisson->SetLineWidth(2);
1014 
1015  return datapoisson;
1016 
1017  }
1018 
1019  //----------------------------------------------------------------------
1020  TGraph* graphAsymmErrorWithBkgScaled(TH1* data, TH1* bkgd, double overallScale){
1021 
1022  for(int bin = 1; bin <= data->GetNbinsX(); bin++) {
1023  if(data->GetBinContent(bin)==0){
1024  data->SetBinContent(bin,0.001);
1025  }
1026  }
1027 
1028  TFeldmanCousins fc(0.6827); // to get correct poisson error bars
1029  TGraphAsymmErrors *gr = new TGraphAsymmErrors(data);
1030 
1031  for(int bin = 1; bin <= data->GetNbinsX(); bin++) {
1032  //default numu: scale by bin 0.1/width
1033  double width = data->GetBinWidth(bin);
1034  double scale = width / overallScale; // this is to reverse scaling
1035 
1036  double x, y;
1037  gr->GetPoint(bin-1, x, y);
1038  double bkg = bkgd->GetBinContent(bin);
1039  y *= scale;
1040  bkg *= scale;
1041 
1042  // since this is data it should be a whole number
1043  // fc Upper/Lower limit is finicky eg if y = 3.01
1044  // the returned value is the same as when y = 4
1045  y = round(y);
1046 
1047  double errup = 0, errdn = 0;
1048  if(y < 30) errup = (fc.CalculateUpperLimit(y,bkg) - y);
1049  else errup = sqrt(y);
1050  if(y < 50) errdn = (y - fc.CalculateLowerLimit(y,bkg));
1051  else errdn = sqrt(y);
1052 
1053  gr->SetPointEYhigh(bin-1, 1/scale * errup);
1054  gr->SetPointEYlow (bin-1, 1/scale * errdn);
1055  }
1056 
1057  gr->SetMarkerStyle(kFullCircle);
1058  gr->SetMarkerColor(1);
1059  gr->SetMarkerSize(1);
1060  gr->SetLineWidth(2);
1061 
1062  return gr;
1063 
1064  }
1065 
1066  //----------------------------------------------------------------------
1067  TGraph* RatioAsymmError(TH1* data, TH1* pred)
1068  {
1069  TGraphAsymmErrors* gRat = GraphWithPoissonErrors(data);
1070  TFeldmanCousins fc(0.6827); // to get correct poisson error bars
1071 
1072  // Hist underflow bin is 0 so index from 1
1073  // Graphs start at 0 though so will offset points
1074  for (int bin = 1; bin <= data->GetNbinsX(); ++bin) {
1075  double x, y;
1076  gRat->GetPoint(bin-1, x, y);
1077  double den = pred->GetBinContent(bin);
1078 
1079  // since this is data it should be a whole number
1080  // fc Upper/Lower limit is finicky eg if y = 3.01
1081  // the returned value is the same as when y = 4
1082  // Also negatives are weird
1083  y = round(y);
1084  if (y <= 0) y = 0.01;
1085 
1086  double errup = 0, errdn = 0;
1087  if (y < 30) errup = (fc.CalculateUpperLimit(y, 0) - y) / den;
1088  else errup = sqrt(y) / den;
1089  if (y < 50) errdn = (y - fc.CalculateLowerLimit(y, 0)) / den;
1090  else errdn = sqrt(y) / den;
1091 
1092  gRat->SetPoint(bin-1, x, y / den);
1093  gRat->SetPointEYhigh(bin-1, errup);
1094  gRat->SetPointEYlow(bin-1, errdn);
1095  }
1096  gRat->SetMarkerStyle(kFullCircle);
1097  gRat->SetMarkerColor(1);
1098  gRat->SetMarkerSize(1);
1099  gRat->SetLineWidth(2);
1100 
1101  return gRat;
1102  }
1103 
1104  //----------------------------------------------------------------------
1105  TGraph* RatioAsymmErrorScaled(TH1* data, TH1* pred, double overallScale)
1106  {
1107  TGraphAsymmErrors* gRat = GraphWithPoissonErrors(data);
1108  TFeldmanCousins fc(0.6827); // to get correct poisson error bars
1109 
1110  // Hist underflow bin is 0 so index from 1
1111  // Graphs start at 0 though so will offset points
1112  for (int bin = 1; bin <= data->GetNbinsX(); ++bin) {
1113  double width = data->GetBinWidth(bin);
1114  double scale = width / overallScale; // this is to reverse scaling
1115 
1116  double x, y;
1117  gRat->GetPoint(bin-1, x, y);
1118  double den = pred->GetBinContent(bin);
1119  y *= scale;
1120  den *= scale;
1121 
1122  // since this is data it should be a whole number
1123  // fc Upper/Lower limit is finicky eg if y = 3.01
1124  // the returned value is the same as when y = 4
1125  // Also negatives are weird
1126  y = round(y);
1127  if (y <= 0) y = 0.01;
1128 
1129  double errup = 0, errdn = 0;
1130  if (y < 30) errup = (fc.CalculateUpperLimit(y, 0) - y) / den;
1131  else errup = sqrt(y) / den;
1132  if (y < 50) errdn = (y - fc.CalculateLowerLimit(y, 0)) / den;
1133  else errdn = sqrt(y) / den;
1134 
1135  gRat->SetPoint(bin-1, x, y / den);
1136  gRat->SetPointEYhigh(bin-1, errup);
1137  gRat->SetPointEYlow (bin-1, errdn);
1138  }
1139  gRat->SetMarkerStyle(kFullCircle);
1140  gRat->SetMarkerColor(1);
1141  gRat->SetMarkerSize(1);
1142  gRat->SetLineWidth(2);
1143 
1144  return gRat;
1145  }
1146 
1147  //----------------------------------------------------------------------
1148  TGraph* RatioAsymmErrorWithBkg(TH1* data, TH1* pred, TH1* bkgd)
1149  {
1150  TGraphAsymmErrors* gRat = GraphWithPoissonErrors(data);
1151  TFeldmanCousins fc(0.6827); // to get correct poisson error bars
1152 
1153  // Hist underflow bin is 0 so index from 1
1154  // Graphs start at 0 though so will offset points
1155  for (int bin = 1; bin <= data->GetNbinsX(); ++bin) {
1156  double x, y;
1157  gRat->GetPoint(bin-1, x, y);
1158  double den = pred->GetBinContent(bin);
1159  double bkg = bkgd->GetBinContent(bin);
1160 
1161  // since this is data it should be a whole number
1162  // fc Upper/Lower limit is finicky eg if y = 3.01
1163  // the returned value is the same as when y = 4
1164  // Also negatives are weird
1165  y = round(y);
1166  if (y <= 0) y = 0.01;
1167 
1168  double errup = 0, errdn = 0;
1169  if (y < 30) errup = (fc.CalculateUpperLimit(y, bkg) - y) / den;
1170  else errup = sqrt(y) / den;
1171  if (y < 50) errdn = (y - fc.CalculateLowerLimit(y, bkg)) / den;
1172  else errdn = sqrt(y) / den;
1173 
1174  gRat->SetPoint(bin-1, x, y / den);
1175  gRat->SetPointEYhigh(bin-1, errup);
1176  gRat->SetPointEYlow (bin-1, errdn);
1177  }
1178  gRat->SetMarkerStyle(kFullCircle);
1179  gRat->SetMarkerColor(1);
1180  gRat->SetMarkerSize(1);
1181  gRat->SetLineWidth(2);
1182 
1183  return gRat;
1184  }
1185 
1186 
1187  //----------------------------------------------------------------------
1188  TGraph* RatioAsymmErrorWithBkgScaled(TH1* data, TH1* pred, TH1* bkgd, double overallScale)
1189  {
1190  TGraphAsymmErrors* gRat = GraphWithPoissonErrors(data);
1191  TFeldmanCousins fc(0.6827); // to get correct poisson error bars
1192 
1193  // Hist underflow bin is 0 so index from 1
1194  // Graphs start at 0 though so will offset points
1195  for (int bin = 1; bin <= data->GetNbinsX(); ++bin) {
1196  double width = data->GetBinWidth(bin);
1197  double scale = width / overallScale; // this is to reverse scaling
1198 
1199  double x, y;
1200  gRat->GetPoint(bin-1, x, y);
1201  double den = pred->GetBinContent(bin);
1202  double bkg = bkgd->GetBinContent(bin);
1203 
1204  y *= scale;
1205  den *= scale;
1206  bkg *= scale;
1207 
1208  // since this is data it should be a whole number
1209  // fc Upper/Lower limit is finicky eg if y = 3.01
1210  // the returned value is the same as when y = 4
1211  // Also negatives are weird
1212  y = round(y);
1213  if (y <= 0) y = 0.01;
1214 
1215  double errup = 0, errdn = 0;
1216  if (y < 30) errup = (fc.CalculateUpperLimit(y, bkg) - y) / den;
1217  else errup = sqrt(y) / den;
1218  if (y < 50) errdn = (y - fc.CalculateLowerLimit(y, bkg)) / den;
1219  else errdn = sqrt(y) / den;
1220 
1221  gRat->SetPoint(bin-1, x, y / den);
1222  gRat->SetPointEYhigh(bin-1, errup);
1223  gRat->SetPointEYlow (bin-1, errdn);
1224  }
1225  gRat->SetMarkerStyle(kFullCircle);
1226  gRat->SetMarkerColor(1);
1227  gRat->SetMarkerSize(1);
1228  gRat->SetLineWidth(2);
1229 
1230  return gRat;
1231  }
1232 
1233  //----------------------------------------------------------------------
1234  TGraph* ShadeBetweenHistograms(TH1* hmin, TH1* hmax)
1235  {
1236  int n = hmin->GetNbinsX();
1237  TGraph* gr = new TGraph(4*n);
1238 
1239  for(int i=1; i<=n; i++)
1240  {
1241  double xdown = hmax->GetBinLowEdge(i);
1242  double xup = hmax->GetBinLowEdge(i+1);
1243  double ydown = hmin->GetBinContent(i);
1244  double yup = hmax->GetBinContent(i);
1245 
1246  gr->SetPoint(2*(i-1), xdown, ydown);
1247  gr->SetPoint(2*(i-1)+1, xup, ydown);
1248  gr->SetPoint(4*n-2*i, xup, yup);
1249  gr->SetPoint(4*n-2*i+1, xdown, yup);
1250  }
1251 
1252  gr->SetFillColor(hmin->GetLineColor() - 9); // In principle this is lighter
1253  gr->SetLineColor(hmin->GetLineColor() - 9); // In case one does Draw("lf")
1254 
1255  return gr;
1256  }
1257  //----------------------------------------------------------------------
1258  TGraphAsymmErrors * ProfileQuantile(const TH2 * hist,
1259  const std::string & axisName,
1260  const std::string & graphName,
1261  const std::pair<double, double> & quantileDivisions)
1262  {
1263  if (hist->GetDimension() != 2)
1264  throw std::runtime_error(Form("Can't profile a histogram with other than 2 dimensions. Yours is a %d-D histogram...", hist->GetDimension()));
1265 
1266  TAxis const* axis = nullptr;
1267  std::function<TH1D*(const char *, Int_t, Int_t, Option_t*)> projectionMethod;
1268  if (axisName == "x" || axisName == "X")
1269  {
1270  axis = hist->GetXaxis();
1271  projectionMethod = std::bind(&TH2::ProjectionY, hist, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4);
1272  }
1273  else if (axisName == "y" || axisName == "Y")
1274  {
1275  axis = hist->GetYaxis();
1276  projectionMethod = std::bind(&TH2::ProjectionX, hist, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4);
1277  }
1278  else
1279  throw std::runtime_error( Form("Can't profile onto unknown axis: '%s'", axisName.c_str()) );
1280 
1281 
1282  std::vector<double> points_x, points_y, errors_x_low, errors_x_high, errors_y_low, errors_y_high;
1283 
1284  double quantiles[2] = {0, 0};
1285  double _quantileDivisions[2] = {quantileDivisions.first, quantileDivisions.second};
1286  for (int binNum = 1; binNum <= axis->GetNbins(); binNum++) // remember, 0 is underflow and Nbins+1 is overflow...
1287  {
1288  std::unique_ptr<TH1D> proj ( projectionMethod(Form("%s_tmp_%d", hist->GetName(), int(std::time(nullptr))), binNum, binNum, "") ); // randomish name, just to be safe
1289 
1290  double y = proj->GetMean();
1291  points_x.push_back(axis->GetBinCenter(binNum));
1292  points_y.push_back(y);
1293 
1294  // for now, make the errors for an empty projection 0
1295  unsigned int n_qs = 0;
1296  if (proj->Integral() == 0)
1297  {
1298  n_qs = 2;
1299  quantiles[0] = quantiles[1] = 0;
1300  }
1301  else
1302  n_qs = proj->GetQuantiles(2, quantiles, _quantileDivisions);
1303  if (n_qs != 2)
1304  throw std::runtime_error( Form("GetQuantiles() didn't compute all the quantiles in HistoTools.ProfileQuantile(). I requested 2 quantiles, but got %d of them...", n_qs) );
1305 
1306  double binWidth = axis->GetBinWidth(binNum);
1307  errors_x_low.push_back(binWidth/2.);
1308  errors_x_high.push_back(binWidth/2.);
1309  errors_y_low.push_back(y-quantiles[0]);
1310  errors_y_high.push_back(quantiles[1]-y);
1311  }
1312 
1313  TGraphAsymmErrors * outGraph = new TGraphAsymmErrors(
1314  points_x.size(),
1315  &points_x[0],
1316  &points_y[0],
1317  &errors_x_low[0],
1318  &errors_x_high[0],
1319  &errors_y_low[0],
1320  &errors_y_high[0]
1321  );
1322  std::string name = (graphName.size()) ? graphName : Form("%s_quantile_errors", hist->GetName());
1323  outGraph->SetName(name.c_str());
1324 
1325  return outGraph;
1326 
1327  }
1328 
1329 
1330  // Draw single BF
1331  void drawBFSingle(double bfSin, double bfDm, Color_t color, Style_t marker, double size = 1.5)
1332  {
1333 
1334  TMarker* manMarker = new TMarker(bfSin, bfDm, marker);
1335  manMarker->SetMarkerSize(size);
1336  manMarker->SetMarkerColor(color);
1337  manMarker->Draw();
1338 
1339  }
1340 
1341  // Truth for FHC; not sure about RHC
1342  void drawBFMirror(double bfSin, double bfDm, Color_t color, Style_t marker, double size = 1.5)
1343  {
1344 
1345  double mirror = bfSin;
1346  TMarker* manMarker = new TMarker(bfSin, bfDm, marker);
1347  manMarker->SetMarkerSize(size);
1348  manMarker->SetMarkerColor(color);
1349  manMarker->Draw();
1350 
1351  if(bfSin > 0.514) mirror = 0.514 - (bfSin - 0.514);
1352  if(bfSin < 0.514) mirror = 0.514 + (0.514 - bfSin );
1353  TMarker* mirrorMarker = new TMarker(mirror, bfDm, marker);
1354  mirrorMarker->SetMarkerSize(size);
1355  mirrorMarker->SetMarkerColor(color);
1356  mirrorMarker->Draw();
1357 
1358  }
1359 
1360 
1361 
1362  //----------------------------------------------------------------------
1363  void MakeHistCanvasReady_Quant(const int quant,
1364  TH1* hist,
1365  double ymax){
1366 
1367  hist->GetXaxis()->SetLabelOffset(0.01);
1368  hist->GetYaxis()->SetLabelOffset(0.01);
1369  hist->GetYaxis()->SetTitleOffset(0.50);
1370  hist->GetXaxis()->SetTitleSize(0.0);
1371  hist->GetYaxis()->SetTitleSize(0.0);
1372  hist->GetYaxis()->SetRangeUser(0.0, ymax);
1373 
1374  if(quant==1 || quant==2) hist->GetXaxis()->SetLabelSize(0.);
1375  if(quant==2 || quant==4) hist->GetYaxis()->SetLabelSize(0.);
1376  // Q3 might need some axis labels removed to prevent interference.
1377  // Rather than assuming there will be trouble we can check the
1378  // optimization that ROOT does when plottng the axes
1379  // ChangeLabel labnum -1 alters the last label
1380  if(quant == 3) {
1381  // initialize the optimized min, max, step size,
1382  // and number of steps of the labels
1383  double optmin = 0., optmax = 0., optwi = 0.;
1384  int optstp = 0;
1385 
1386  double xmin = hist->GetXaxis()->GetXmin();
1387  double xmax = hist->GetXaxis()->GetXmax();
1388  int xdiv = hist->GetXaxis()->GetNdivisions() % 100;
1389  THLimitsFinder::Optimize(xmin, xmax, xdiv, optmin, optmax, optstp, optwi);
1390 
1391  if ( (xmax - optmax) < optwi/2 )
1392  hist->GetXaxis()->ChangeLabel( -1, -1, 0);
1393 
1394  optmin = 0.; optmax = 0.; optwi = 0.; optstp = 0;
1395  double ymin = 0.;
1396  int ydiv = hist->GetYaxis()->GetNdivisions() % 100;
1397  THLimitsFinder::Optimize(ymin, ymax, ydiv, optmin, optmax, optstp, optwi);
1398 
1399  hist->GetYaxis()->SetRangeUser(ymin, ymax);
1400  if ( (ymax - optmax) < optwi/2 )
1401  hist->GetYaxis()->ChangeLabel( -1, -1, 0);
1402  }
1403  }
1404 
1405  //----------------------------------------------------------------------
1406  void PimpHist(TH1* hist, Style_t linestyle, Color_t linecolor, int linewidth, Style_t markerstyle, Color_t markercolor, double markersize){
1407  hist->SetLineColor(linecolor);
1408  hist->SetLineStyle(linestyle);
1409  hist->SetLineWidth(linewidth);
1410  hist->SetMarkerColor(markercolor);
1411  hist->SetMarkerStyle(markerstyle);
1412  hist->SetMarkerSize(markersize);
1413  }
1414 
1415 
1416  //----------------------------------------------------------------------
1417  void SplitCanvas(double ysplit, TPad*& p1, TPad*& p2)
1418  {
1419  if(!gPad) new TCanvas;
1420  p1 = new TPad("", "", 0, 0, 1, 1);
1421  p2 = new TPad("", "", 0, 0, 1, 1);
1422 
1423  p1->SetBottomMargin(ysplit);
1424  p2->SetTopMargin(1-ysplit);
1425 
1426  // Draw p1 second since it's often the more important one, that the user
1427  // would prefer to be able to interact with.
1428  for(TPad* p: {p2, p1}){
1429  p->SetFillStyle(0);
1430  p->Draw();
1431  }
1432  }
1433 
1434 
1435  //----------------------------------------------------------------------
1436  // split canvas in 4 for numu quartiles
1437  void SplitCanvasQuant(TCanvas *& canvas, TPad *& pad1, TPad *& pad2, TPad *& pad3, TPad *& pad4){
1438 
1439  canvas = new TCanvas(UniqueName().c_str(),"",960,800);
1440 
1441  pad3 = new TPad(UniqueName().c_str(),"pad3",0,0,1,1);
1442  pad3->Draw();
1443  pad3->SetTopMargin(0.5);
1444  pad3->SetBottomMargin(0.1);
1445  pad3->SetLeftMargin(0.1);
1446  pad3->SetRightMargin(0.49);
1447  pad3->SetFillStyle(0);
1448  canvas->cd();
1449 
1450  pad4 = new TPad(UniqueName().c_str(),"pad4",0,0,1,1);
1451  pad4->Draw();
1452  pad4->SetTopMargin(0.5);
1453  pad4->SetBottomMargin(0.1);
1454  pad4->SetLeftMargin(0.51);
1455  pad4->SetRightMargin(0.08);
1456  pad4->SetFillStyle(0);
1457  canvas->cd();
1458 
1459  pad1 = new TPad(UniqueName().c_str(),"pad1",0,0,1,1);// x1 y1 x2 y2
1460  pad1->Draw();
1461  pad1->SetTopMargin(0.1);
1462  pad1->SetBottomMargin(0.5);
1463  pad1->SetLeftMargin(0.1);
1464  pad1->SetRightMargin(0.49);
1465  pad1->SetFillStyle(0);
1466  canvas->cd();
1467 
1468  pad2 = new TPad(UniqueName().c_str(),"pad2",0,0,1,1);
1469  pad2->Draw();
1470  pad2->SetTopMargin(0.1);
1471  pad2->SetBottomMargin(0.5);
1472  pad2->SetLeftMargin(0.51);
1473  pad2->SetRightMargin(0.08);
1474  pad2->SetFillStyle(0);
1475  canvas->cd();
1476 
1477  }
1478 
1479 
1480  //----------------------------------------------------------------------
1481  void CenterTitles ( TH1 * histo)
1482  {
1483  histo->GetXaxis()->CenterTitle();
1484  histo->GetYaxis()->CenterTitle();
1485  histo->GetZaxis()->CenterTitle();
1486  }
1487 
1488 
1489  //----------------------------------------------------------------------
1490  bool SortSystsName(const ISyst* systA, const ISyst* systB){
1491  return ( systA->ShortName() < systB->ShortName() );
1492  }
1493 
1494 
1495  TH1* PlotSystShifts(const SystShifts & shifts, bool sortName){
1496 
1497  auto systs = shifts.ActiveSysts();
1498  int nsysts = systs.size();
1499 
1500  // sort systematics by name
1501  if(sortName) std::sort(systs.begin(), systs.end(), SortSystsName);
1502 
1503  auto h = new TH1D ("sh",";; Pull (N#sigma)",
1504  nsysts, -0.5, nsysts + 0.5);
1505 
1506  for (int systIdx = 0; systIdx < nsysts; ++systIdx){
1507  double shiftThis = shifts.GetShift(systs[systIdx]);
1508  h->SetBinContent( systIdx + 1, shiftThis);
1509 
1510  TString tempName = systs[systIdx]->ShortName();
1511  if( tempName.Contains("Neutron") ) h->GetXaxis()->SetBinLabel( systIdx + 1, "NeutronSyst2018");
1512  else h->GetXaxis()->SetBinLabel( systIdx + 1, systs[systIdx]->ShortName().c_str());
1513  }
1514 
1515  h->SetMarkerStyle(kFullCircle);
1516  CenterTitles(h);
1517  h->SetTitleOffset(3.0);
1518  h->GetXaxis()->SetNdivisions(nsysts,kFALSE);
1519  h->GetXaxis()->SetLabelSize(0.065);
1520  h->GetXaxis()->LabelsOption("v");
1521  h->GetYaxis()->SetRangeUser(-1,1);
1522  auto c1 = new TCanvas ("c1","c1",700,350);
1523  c1->SetBottomMargin(0.5);
1524  h->Draw("lp hist");
1525  return h;
1526  }
1527 
1528  //----------------------------------------------------------------------
1529  TGraph* JoinGraphs(TGraph* a, TGraph* b, int fillcol)
1530  {
1531  TGraph* ret = (TGraph*)a->Clone(UniqueName().c_str());
1532  for(int i = 0; i < b->GetN(); ++i){
1533  double x, y;
1534  b->GetPoint(i, x, y);
1535  ret->SetPoint(ret->GetN(), x, y);
1536  }
1537 
1538  ret->SetFillColor(fillcol);
1539 
1540  return ret;
1541  }
1542 
1543  //----------------------------------------------------------------------
1544  TGraph* ExtendGraphToTop(TGraph* g, int col, double xmin, double xmax, double y)
1545  {
1546  g->SetPoint(g->GetN(), xmax, y);
1547  g->SetPoint(g->GetN(), xmin, y);
1548 
1549  g->SetFillColor(col);
1550 
1551  return g;
1552  }
1553 
1554  //----------------------------------------------------------------------
1555  TPaveText* DrawBeamLabel(bool isFHC)
1556  {
1557 
1558  TPaveText *pText = new TPaveText (0.1, 0.90, 0.16, 0.94,"NDC");
1559  pText->SetFillStyle(0);
1560  pText->SetLineColor(0);
1561  pText->SetLineWidth(0);
1562  pText->SetBorderSize(1);
1563  pText->SetTextColor(kGray+2);
1564  if (isFHC) pText->AddText("#nu-beam");
1565  else pText->AddText("#bar{#nu}-beam");
1566  pText->SetTextSize(2/40.);
1567  pText->SetTextAlign(11);
1568  pText->Draw();
1569 
1570  return pText;
1571  }
1572 
1573  //----------------------------------------------------------------------
1574  TPaveText* DrawQuantLabel(int quant)
1575  {
1576  TPaveText *ret = 0;
1577  // Draw nothing for q = 0
1578  if (quant == 0) {
1579  // Draw nothing for q = 0 but it's still valid
1580  }
1581  else if (quant == 1) {
1582  ret = new TPaveText(0.37, 0.75, 0.47, 0.90, "NDC");
1583  ret->SetBorderSize(0);
1584  ret->SetFillStyle(0);
1585  TText *text = ret->AddText(TString::Format("#splitline{#bf{Quartile %d}}{best resolution}", quant));
1586  text ->SetTextAlign(22);
1587  text ->SetTextSize(0.025);
1588  ret->Draw();
1589  }
1590  else if (quant == 2) {
1591  ret = new TPaveText(0.75, 0.75, 0.85, 0.90, "NDC");
1592  ret->SetBorderSize(0);
1593  ret->SetFillStyle(0);
1594  TText *text = ret->AddText(TString::Format("#bf{Quartile %d}", quant));
1595  text ->SetTextAlign(22);
1596  text ->SetTextSize(0.025);
1597  ret->Draw();
1598  }
1599  else if (quant == 3) {
1600  ret = new TPaveText(0.37, 0.35, 0.47, 0.50, "NDC");
1601  ret->SetBorderSize(0);
1602  ret->SetFillStyle(0);
1603  TText *text = ret->AddText(TString::Format("#bf{Quartile %d}", quant));
1604  text ->SetTextAlign(22);
1605  text ->SetTextSize(0.025);
1606  ret->Draw();
1607  }
1608  else if (quant == 4) {
1609  ret = new TPaveText(0.75, 0.35, 0.85, 0.50, "NDC");
1610  ret->SetBorderSize(0);
1611  ret->SetFillStyle(0);
1612  TText *text = ret->AddText(TString::Format("#splitline{#bf{Quartile %d}}{worst resolution}", quant));
1613  text ->SetTextAlign(22);
1614  text ->SetTextSize(0.025);
1615  ret->Draw();
1616  }
1617  else {
1618  std::cout << "\n***********************************************\n"
1619  << "\n No valid quantile chosen for DrawQuantLabel \n"
1620  << "\n***********************************************\n"
1621  << std::endl;
1622  }
1623 
1624  return ret;
1625  }
1626 
1627 } // namespace
T max(const caf::Proxy< T > &a, T b)
void PlotWithAreaSystErrorBand(IPrediction *pred, const std::vector< const ISyst * > &systs, osc::IOscCalc *calc, double pot, int col, int errCol, float headroom, bool newaxis, EBinType bintype)
Plot prediction with +/-1sigma shape-only error band.
Definition: Plots.cxx:417
const XML_Char * name
Definition: expat.h:151
std::vector< double > quantiles(TH1D *h)
Definition: absCal.cxx:528
TSpline3 lo("lo", xlo, ylo, 12,"0")
void DataMCRatio(const Spectrum &data, const IPrediction *mc, osc::IOscCalc *calc, double miny, double maxy)
Plot data/MC ratio for the given spectrum. Normalize MC to Data by POT.
Definition: Plots.cxx:167
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::map< std::string, double > xmax
virtual Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
set< int >::iterator it
TH1 * PlotSystShifts(const SystShifts &shifts, bool sortName)
Definition: Plots.cxx:1495
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
TH1D * ToTH1(Color_t col=kBlack, Style_t style=kSolid) const
Definition: Ratio.cxx:68
void MakeHistCanvasReady_Quant(const int quant, TH1 *hist, double ymax)
Definition: Plots.cxx:1363
Eigen::ArrayXd ProjectionX(const Eigen::MatrixXd &mat)
Helper for Unweighted.
Float_t y1[n_points_granero]
Definition: compare.C:5
bool IsNominal() const
Definition: SystShifts.h:43
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:209
virtual const std::string & ShortName() const final
The name printed out to the screen.
Definition: ISyst.h:27
virtual Spectrum Predict(osc::IOscCalc *calc) const =0
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
void OverridePOT(double newpot)
DO NOT USE UNLESS YOU ARE 110% CERTAIN THERE ISN&#39;T A BETTER WAY!
Definition: Spectrum.h:237
double maxy
Float_t x1[n_points_granero]
Definition: compare.C:5
Float_t den
Definition: plot.C:36
(&#39;beam &#39;)
Definition: IPrediction.h:15
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
TH1 * ratio(TH1 *h1, TH1 *h2)
General interface to oscillation calculators.
Definition: FwdDeclare.h:14
virtual Spectrum PredictComponentSyst(osc::IOscCalc *calc, const SystShifts &syst, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const
TPaveText * DrawQuantLabel(int quant)
Definition: Plots.cxx:1574
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
void CenterTitles(TH1 *histo)
Definition: Plots.cxx:1481
TGraphAsymmErrors * GraphWithPoissonErrors(const TH1 *h, bool noErrorsXaxis, bool drawEmptyBins)
Calculate statistical errors appropriate for small Poisson numbers.
Definition: Plots.cxx:910
void SplitCanvas(double ysplit, TPad *&p1, TPad *&p2)
Split the current pad into two vertically stacked pieces, suitable for ratio plots and so on...
Definition: Plots.cxx:1417
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:333
const Color_t kTotalMCErrorBandColor
Definition: Style.h:17
static SystShifts Nominal()
Definition: SystShifts.h:34
fvar< T > round(const fvar< T > &x)
Definition: round.hpp:23
osc::OscCalcDumb calc
::xsd::cxx::tree::time< char, simple_type > time
Definition: Database.h:194
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
const Color_t kNumuBackgroundColor
Definition: Style.h:30
void GetSystBands(osc::IOscCalc *calc, const IPrediction *pred, std::vector< const ISyst * > systs, std::vector< TH1 * > &hUps, std::vector< TH1 * > &hDns, double pot, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign)
Definition: Plots.cxx:244
TGraph * RatioAsymmError(TH1 *data, TH1 *pred)
Definition: Plots.cxx:1067
TGraphAsymmErrors * ProfileQuantile(const TH2 *hist, const std::string &axisName, const std::string &graphName, const std::pair< double, double > &quantileDivisions)
Calculate profile with error bars corresponding to specified quantiles of a 2D distribution (by defau...
Definition: Plots.cxx:1258
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
Float_t Y
Definition: plot.C:38
const char * label
TGraphAsymmErrors * PlotWithSystErrorBand(IPrediction *pred, const std::vector< const ISyst * > &systs, osc::IOscCalc *calc, double pot, int col, int errCol, float headroom, bool newaxis, EBinType bintype)
Plot prediction with +/-1sigma error band.
Definition: Plots.cxx:304
TGraph * gr2
Definition: compare.C:43
void PimpHist(TH1 *hist, Style_t linestyle, Color_t linecolor, int linewidth, Style_t markerstyle, Color_t markercolor, double markersize)
Pimp histogram once and for all.
Definition: Plots.cxx:1406
const XML_Char const XML_Char * data
Definition: expat.h:268
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
Double_t ymax
Definition: plot.C:25
EBinType
Definition: Utilities.h:42
void drawBFSingle(double bfSin, double bfDm, Color_t color, Style_t marker, double size=1.5)
Draw a single BF point.
Definition: Plots.cxx:1331
const XML_Char * s
Definition: expat.h:262
const int nbins
Definition: cellShifts.C:15
Double_t scale
Definition: plot.C:25
Charged-current interactions.
Definition: IPrediction.h:39
TSpline3 hi("hi", xhi, yhi, 18,"0")
double dy[NP][NC]
double dx[NP][NC]
void PlotWithSystErrorBandTwoPreds(const Spectrum &nominal, const std::vector< Spectrum > &upShifts, const std::vector< Spectrum > &downShifts, const Spectrum &alternative, const std::vector< Spectrum > &upShiftsAlt, const std::vector< Spectrum > &downShiftsAlt, double pot, int col1, int col2, float headroom, bool newaxis, EBinType bintype)
Plot two different predictions with +/-1sigma shape-only error bands.
Definition: Plots.cxx:468
TH1 * DataMCComparisonAreaNormalized(const Spectrum &data, const Spectrum &mc)
Definition: Plots.cxx:74
TGraph * graphAsymmErrorScaled(TH1 *histScaled, TH1 *hist, double overallScale)
Definition: Plots.cxx:978
TPaveText * DrawBeamLabel(bool isFHC)
Put the standardized beam label in the left corner of the active canvas.
Definition: Plots.cxx:1555
Int_t col[ntarg]
Definition: Style.C:29
const double a
const Color_t kBeamNueBackgroundColor
Definition: Style.h:24
TGraph * ShadeBetweenHistograms(TH1 *hmin, TH1 *hmax)
Definition: Plots.cxx:1234
TGraphAsymmErrors * GraphWithPoissonErrors2(const TH1 *h, const TH1 *h2, bool noErrorsXaxis, bool drawEmptyBins)
Same as above but use a reference histogram to determine which empty bins to draw.
Definition: Plots.cxx:942
#define pot
T GetShift(const ISyst *syst) const
Float_t d
Definition: plot.C:236
virtual Spectrum PredictSyst(osc::IOscCalc *calc, const SystShifts &syst) const
Definition: IPrediction.cxx:72
TGraph * ExtendGraphToTop(TGraph *g, int col, double xmin, double xmax, double y)
Definition: Plots.cxx:1544
const double j
Definition: BetheBloch.cxx:29
TH2D * histo
float bin[41]
Definition: plottest35.C:14
double POT() const
Definition: Spectrum.h:231
static bool isFHC
Represent the ratio between two spectra.
Definition: Ratio.h:8
TGraph * JoinGraphs(TGraph *a, TGraph *b, int fillcol)
Join graphs and set a fill color. Useful for contours.
Definition: Plots.cxx:1529
OStream cout
Definition: OStream.cxx:6
TH1F * h2
Definition: plot.C:45
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
double PointDistanceToBox(double x, double y, double x0, double y0, double x1, double y1)
Helper for AutoPlaceLegend.
Definition: Plots.cxx:688
TGraph * RatioAsymmErrorWithBkg(TH1 *data, TH1 *pred, TH1 *bkgd)
Definition: Plots.cxx:1148
void CountingExperimentErrorBarChart(const std::map< std::string, double > &systbars, double statErr, bool bkgdOrSig, bool shortchart)
Make a simple plot of relative size of different errors.
Definition: Plots.cxx:896
void DataMCAreaNormalizedRatio(const Spectrum &data, const IPrediction *mc, osc::IOscCalc *calc, double miny, double maxy)
Plot data/MC ratio for the given spectrum. Normalize MC to Data by area.
Definition: Plots.cxx:196
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
TGraph * graphAsymmErrorWithBkgScaled(TH1 *data, TH1 *bkgd, double overallScale)
Definition: Plots.cxx:1020
TH1 * DataMCComparisonComponents(const Spectrum &data, const IPrediction *mc, osc::IOscCalc *calc)
Plot MC broken down into flavour components, overlayed with data.
Definition: Plots.cxx:114
Eigen::ArrayXd ProjectionY(const Eigen::MatrixXd &mat)
Helper for WeightingVariable.
TGraphAsymmErrors * PlotWithSystErrorBand_Quant(const int quant, IPrediction *pred, const std::vector< const ISyst * > &systs, osc::IOscCalc *calc, double pot, int col, int errCol, float maxy, bool newaxis)
Definition: Plots.cxx:541
void GetBFSystBands(osc::IOscCalc *calc, const IPrediction *pred, std::vector< const ISyst * > systs, const SystShifts &bfshifts, std::vector< TH1 * > &hUps, std::vector< TH1 * > &hDns, double pot, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign)
Definition: Plots.cxx:271
Float_t proj
Definition: plot.C:35
TGraph * RatioAsymmErrorWithBkgScaled(TH1 *data, TH1 *pred, TH1 *bkgd, double overallScale)
Definition: Plots.cxx:1188
const hit & b
Definition: hits.cxx:21
Neutral-current interactions.
Definition: IPrediction.h:40
void ErrorBarChart(const std::map< std::string, std::pair< double, double >> &systErrors, const std::pair< double, double > &statErr, const std::string &label)
Make a simple plot of relative size of different errors.
Definition: Plots.cxx:798
THStack * ToTHStack(const std::vector< std::pair< const Spectrum &, Color_t >> &s, double pot)
Can call like ToTHStack({{h1, kRed}, {h2, kBlue}}, pot)
Definition: Plots.cxx:674
bool SortSystsName(const ISyst *systA, const ISyst *systB)
Definition: Plots.cxx:1490
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
c1
Definition: demo5.py:24
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
Double_t ymin
Definition: plot.C:24
TPad * pad3
Definition: analysis.C:13
std::vector< const ISyst * > ActiveSysts() const
Definition: SystShifts.cxx:218
string tempName
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
const Color_t kTotalMCColor
Definition: Style.h:16
void SplitCanvasQuant(TCanvas *&canvas, TPad *&pad1, TPad *&pad2, TPad *&pad3, TPad *&pad4)
Definition: Plots.cxx:1437
TGraph * RatioAsymmErrorScaled(TH1 *data, TH1 *pred, double overallScale)
Definition: Plots.cxx:1105
auto one()
Definition: PMNS.cxx:49
TPad * pad2
Definition: analysis.C:13
TH1 * DataMCComparison(const Spectrum &data, const Spectrum &mc, EBinType bintype)
Definition: Plots.cxx:35
T min(const caf::Proxy< T > &a, T b)
All neutrinos, any flavor.
Definition: IPrediction.h:26
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
Float_t X
Definition: plot.C:38
const Color_t kNCBackgroundColor
Definition: Style.h:22
Float_t w
Definition: plot.C:20
TCanvas * RatioPlot(std::vector< TH1 * > topHistos, std::vector< TString > topOption, std::vector< TH1 * > bottomHistos, std::vector< TString > bottomOption, TString pidtag, bool pidaxis=false)
auto zero()
Definition: PMNS.cxx:47
void next()
Definition: show_event.C:84
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:78
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:30
void drawBFMirror(double bfSin, double bfDm, Color_t color, Style_t marker, double size=1.5)
Draw best fit at both octants. Truth for neutrinos, not sure about antineutrinos. ...
Definition: Plots.cxx:1342
def sign(x)
Definition: canMan.py:197
TPad * pad1
Definition: analysis.C:13