Plots.cxx
Go to the documentation of this file.
3 
6 #include "CAFAna/Core/Ratio.h"
7 #include "CAFAna/Core/Spectrum.h"
9 
10 #include "Utilities/func/MathUtil.h"
11 
12 #include "TCanvas.h"
13 #include "TFeldmanCousins.h"
14 #include "TGraph.h"
15 #include "TGraphAsymmErrors.h"
16 #include "TH1.h"
17 #include "TH2.h"
18 #include "THStack.h"
19 #include "TLegend.h"
20 #include "TLine.h"
21 #include "TList.h"
22 #include "TMarker.h"
23 
24 #include <algorithm>
25 #include <iostream>
26 #include <cmath>
27 #include <ctime>
28 #include <functional>
29 #include <memory>
30 
31 namespace ana
32 {
33  //----------------------------------------------------------------------
34  TH1* DataMCComparison(const Spectrum& data, const Spectrum& mc, EBinType bintype)
35  {
36  TH1* ret = 0;
37 
38  TH1* hMC = mc.ToTH1(data.POT(), kTotalMCColor, kSolid, kPOT, bintype);
39 
40  TH1* hData = data.ToTH1(data.POT(), kBlack, kSolid, kPOT, bintype);
41  hData->Sumw2();
42  hData->SetMarkerStyle(kFullCircle);
43 
44  // Need to draw the biggest thing first to get correct axis limits
45  if(hMC->GetMaximum() > hData->GetMaximum()+sqrt(hData->GetMaximum())){
46  ret = hMC;
47  hMC->Draw("hist");
48  hData->Draw("ep same");
49  }
50  else{
51  ret = hData;
52  hData->Draw("ep");
53  hMC->Draw("hist same");
54  hData->Draw("ep same"); // Data always goes on top
55  }
56 
57  gPad->Update();
58 
59  return ret;
60  }
61 
62  //----------------------------------------------------------------------
64  const IPrediction* mc,
66  const SystShifts & shifts,
67  EBinType bintype)
68  {
69  return DataMCComparison(data, shifts.IsNominal() ? mc->Predict(calc) : mc->PredictSyst(calc, shifts), bintype);
70  }
71 
72  //----------------------------------------------------------------------
74  {
75  TH1* ret = 0;
76 
77  TH1* hMC = mc.ToTH1(mc.POT());
78  hMC->SetLineColor(kTotalMCColor);
79 
80  TH1* hData = data.ToTH1(data.POT());
81  hData->Sumw2();
82  hData->SetMarkerStyle(kFullCircle);
83 
84  hMC->Scale(hData->Integral()/hMC->Integral());
85 
86  // Need to draw the biggest thing first to get correct axis limits
87  if(hMC->GetMaximum() > hData->GetMaximum()+sqrt(hData->GetMaximum())){
88  ret = hMC;
89  hMC->Draw("hist");
90  hData->Draw("ep same");
91  }
92  else{
93  ret = hData;
94  hData->Draw("ep");
95  hMC->Draw("hist same");
96  hData->Draw("ep same"); // Data always goes on top
97  }
98 
99  gPad->Update();
100 
101  return ret;
102  }
103 
104  //----------------------------------------------------------------------
106  const IPrediction* mc,
108  {
109  return DataMCComparisonAreaNormalized(data, mc->Predict(calc));
110  }
111 
112  //----------------------------------------------------------------------
114  const IPrediction* mc,
116  {
117  TH1* ret = 0;
118 
119  TH1* hMC = mc->Predict(calc).ToTH1(data.POT());
120  hMC->SetLineColor(kTotalMCColor);
121 
122  TH1* hNC = mc->PredictComponent(calc,
124  Current::kNC,
125  Sign::kBoth).ToTH1(data.POT());
126  hNC->SetLineColor(kNCBackgroundColor);
127 
128  TH1* hCC = mc->PredictComponent(calc,
130  Current::kCC,
131  Sign::kBoth).ToTH1(data.POT());
132  hCC->SetLineColor(kNumuBackgroundColor);
133 
134  TH1* hBeam = mc->PredictComponent(calc,
136  Current::kCC,
137  Sign::kBoth).ToTH1(data.POT());
138  hBeam->SetLineColor(kBeamNueBackgroundColor);
139 
140  TH1* hData = data.ToTH1(data.POT());
141  hData->SetMarkerStyle(kFullCircle);
142 
143  // Need to draw the biggest thing first to get correct axis limits
144  if(hMC->GetMaximum() > hData->GetMaximum()+sqrt(hData->GetMaximum())){
145  ret = hMC;
146  hMC->Draw("hist");
147  hData->Draw("ep same");
148  }
149  else{
150  ret = hData;
151  hData->Draw("ep");
152  hMC->Draw("hist same");
153  hData->Draw("ep same"); // Data always goes on top
154  }
155 
156  hNC->Draw("hist same");
157  hCC->Draw("hist same");
158  hBeam->Draw("hist same");
159 
160  gPad->Update();
161 
162  return ret;
163  }
164 
165  //----------------------------------------------------------------------
166  void DataMCRatio(const Spectrum& data,
167  const IPrediction* mc,
169  double miny, double maxy)
170  {
171  DataMCRatio(data, mc->Predict(calc), miny, maxy);
172  }
173 
174  //----------------------------------------------------------------------
175  void DataMCRatio(const Spectrum& data,
176  const Spectrum& mc,
177  double miny, double maxy)
178  {
179  Ratio ratio(data, mc);
180  TH1* h = ratio.ToTH1();
181  h->GetYaxis()->SetTitle("Data / MC");
182  h->GetYaxis()->SetRangeUser(miny, maxy);
183  h->Draw();
184 
185  TLine* one = new TLine(h->GetXaxis()->GetXmin(), 1,
186  h->GetXaxis()->GetXmax(), 1);
187  one->SetLineWidth(2);
188  one->SetLineStyle(7);
189  one->Draw();
190 
191  gPad->Update();
192  }
193 
194  //----------------------------------------------------------------------
196  const IPrediction* mc,
198  double miny, double maxy)
199  {
200  DataMCAreaNormalizedRatio(data, mc->Predict(calc), miny, maxy);
201  }
202 
203  //----------------------------------------------------------------------
205  const Spectrum& mc,
206  double miny, double maxy)
207  {
208  Spectrum mcScaled = mc;
209  mcScaled.OverridePOT(mcScaled.POT() * mc.Integral(1) / data.Integral(1));
210 
211  DataMCRatio(data, mcScaled, miny, maxy);
212  }
213 
214  //----------------------------------------------------------------------
215  void RatioPlot(const Spectrum& data,
216  const Spectrum& expected,
217  const Spectrum& fit,
218  double miny, double maxy)
219  {
220  Ratio fitRatio(fit, expected);
221  Ratio dataRatio(data, expected);
222 
223  TH1* h = fitRatio.ToTH1();
224  h->GetYaxis()->SetTitle("Ratio to expectation");
225  h->GetYaxis()->SetRangeUser(miny, maxy);
226  h->SetLineColor(kTotalMCColor);
227  h->Draw("][");
228 
229  h = dataRatio.ToTH1();
230  h->SetMarkerStyle(kFullCircle);
231  h->Draw("ep same");
232 
233  TLine* one = new TLine(h->GetXaxis()->GetXmin(), 1,
234  h->GetXaxis()->GetXmax(), 1);
235  one->SetLineWidth(2);
236  one->SetLineStyle(7);
237  one->Draw();
238 
239  gPad->Update();
240  }
241 
242 
243  //----------------------------------------------------------------------
245  const std::vector<const ISyst*>& systs,
247  double pot,
248  int col, int errCol, float headroom,
249  bool newaxis,
250  EBinType bintype)
251  {
252 
253  Spectrum nom = pred->Predict(calc);
254 
255  std::vector<Spectrum> ups, dns;
256 
257  for(const ISyst* syst: systs){
258  SystShifts shifts;
259  shifts.SetShift(syst, +1);
260  ups.push_back(pred->PredictSyst(calc, shifts));
261  shifts.SetShift(syst, -1);
262  dns.push_back(pred->PredictSyst(calc, shifts));
263  }
264 
265  return PlotWithSystErrorBand(nom, ups, dns, pot, col, errCol, headroom, newaxis, bintype);
266 
267 
268  }
269 
270  //----------------------------------------------------------------------
271 
272  TGraphAsymmErrors* PlotWithSystErrorBand(const Spectrum& nominal,
273  const std::vector<Spectrum>& upShifts,
274  const std::vector<Spectrum>& downShifts,
275  double pot, int col, int errCol,
276  float headroom, bool newaxis,
277  EBinType bintype)
278  {
279 
280  TH1* nom = nominal.ToTH1(pot, kBlack, kSolid, kPOT, bintype);
281 
282  std::vector<TH1*> ups, dns;
283 
284  for(const auto& upShift:upShifts) ups.push_back(upShift.ToTH1(pot, kBlack, kSolid, kPOT, bintype));
285  for(const auto& downShift:downShifts) dns.push_back(downShift.ToTH1(pot, kBlack, kSolid, kPOT, bintype));
286 
287  return PlotWithSystErrorBand(nom, ups, dns, col, errCol, headroom, newaxis);
288  }
289 
290  //----------------------------------------------------------------------
291 
292  TGraphAsymmErrors* PlotWithSystErrorBand(TH1*& nom,
293  std::vector<TH1*>& ups,
294  std::vector<TH1*>& dns,
295  int col, int errCol,
296  float headroom, bool newaxis)
297  {
298  if(col == -1){
299  col = kTotalMCColor;
300  errCol = kTotalMCErrorBandColor;
301  }
302  else if(errCol == -1) errCol = col-7; // hopefully a lighter version
303 
304  nom->SetLineColor(col);
305  nom->GetXaxis()->CenterTitle();
306  nom->GetYaxis()->CenterTitle();
307  if(newaxis) nom->Draw("hist ]["); // Set the axes up
308 
309  double yMax = nom->GetBinContent(nom->GetMaximumBin());
310 
311  TGraphAsymmErrors* g = new TGraphAsymmErrors;
312 
313  for(int binIdx = 0; binIdx < nom->GetNbinsX()+2; ++binIdx){
314  const double y = nom->GetBinContent(binIdx);
315  g->SetPoint(binIdx, nom->GetXaxis()->GetBinCenter(binIdx), y);
316 
317  const double w = nom->GetXaxis()->GetBinWidth(binIdx);
318 
319  double errUp = 0;
320  double errDn = 0;
321 
322  for(unsigned int systIdx = 0; systIdx < ups.size(); ++systIdx){
323  double hi = ups[systIdx]->GetBinContent(binIdx)-y;
324  double lo = dns[systIdx]->GetBinContent(binIdx)-y;
325 
326  // If both shifts are in the same direction use the larger
327  double min = std::min(hi,lo);
328  double max = std::max(hi,lo);
329  if(max < 0) max=0;
330  if(min > 0) min=0;
331 
332  errUp += max*max;
333  errDn += min*min;
334  } // end for systIdx
335 
336  g->SetPointError(binIdx, w/2, w/2, sqrt(errDn), sqrt(errUp));
337  } // end for i
338 
339  g->SetFillColor(errCol);
340  g->Draw("e2 same");
341  g->GetYaxis()->SetRangeUser(0.00001, headroom*yMax);
342  nom->GetYaxis()->SetRangeUser(0.0001, headroom*yMax);
343 
344  nom->Draw("hist ][ same");
345 
346  for(TH1* up: ups) delete up;
347  for(TH1* dn: dns) delete dn;
348  return g;
349  }
350 
351  //----------------------------------------------------------------------
353  const std::vector<const ISyst*>& systs,
355  double pot,
356  int col, int errCol, float headroom,
357  bool newaxis, EBinType bintype)
358  {
359 
360  Spectrum nom = pred->Predict(calc);
361 
362  std::vector<Spectrum> ups, dns;
363 
364  for(const ISyst* syst: systs){
365  SystShifts shifts;
366  shifts.SetShift(syst, +1);
367  ups.push_back(pred->PredictSyst(calc, shifts));
368  shifts.SetShift(syst, -1);
369  dns.push_back(pred->PredictSyst(calc, shifts));
370  }
371 
372  PlotWithAreaSystErrorBand(nom, ups, dns, pot, col, errCol, headroom, newaxis, bintype);
373 
374 
375  }
376 
377  //----------------------------------------------------------------------
378  void PlotWithAreaSystErrorBand(const Spectrum& nominal,
379  std::vector<Spectrum> upShifts,
380  std::vector<Spectrum> downShifts,
381  double pot, int col, int errCol,
382  float headroom, bool newaxis,
383  EBinType bintype)
384  {
385  double nomIntegral = nominal.Integral(pot);
386  if (nomIntegral > 0)
387  {
388  // ugh so much boilerplate
389  for (auto & collection : std::vector<std::reference_wrapper<std::vector<Spectrum>>>{ upShifts, downShifts })
390  {
391  for (Spectrum & shiftSpec : collection.get())
392  {
393  // this is a dirty trick, but these Spectra are going to be thrown out once this method exits anyway
394  shiftSpec.OverridePOT(shiftSpec.POT() * shiftSpec.Integral(pot) / nomIntegral);
395  }
396  }
397  }
398 
399  PlotWithSystErrorBand(nominal, upShifts, downShifts, pot, col, errCol, headroom, newaxis, bintype);
400  }
401 
402  //----------------------------------------------------------------------
404  const std::vector<Spectrum>& upShifts,
405  const std::vector<Spectrum>& downShifts,
406  const Spectrum& alternative,
407  const std::vector<Spectrum>& upShiftsAlt,
408  const std::vector<Spectrum>& downShiftsAlt,
409  double pot, int col1, int col2,
410  float headroom, bool newaxis,
411  EBinType bintype)
412  {
413  TH1* nom = nominal.ToTH1(pot, kBlack, kSolid, kPOT, bintype);
414  TH1* alt = alternative.ToTH1(pot, kBlack, kSolid, kPOT, bintype);
415  std::vector<TH1*> ups, dns;
416  std::vector<TH1*> upsA, dnsA;
417  for(const auto& upShift:upShifts) ups.push_back(upShift.ToTH1(pot, kBlack, kSolid, kPOT, bintype));
418  for(const auto& downShift:downShifts) dns.push_back(downShift.ToTH1(pot, kBlack, kSolid, kPOT, bintype));
419  for(const auto& upShiftAlt:upShiftsAlt) upsA.push_back(upShiftAlt.ToTH1(pot, kBlack, kSolid, kPOT, bintype));
420  for(const auto& downShiftAlt:downShiftsAlt) dnsA.push_back(downShiftAlt.ToTH1(pot, kBlack, kSolid, kPOT, bintype));
421  nom->SetLineColor(col1);
422  nom->GetXaxis()->CenterTitle();
423  nom->GetYaxis()->CenterTitle();
424  if(newaxis) nom->Draw("hist"); // Set the axes up
425  alt->SetLineColor(col2);
426  double yMax = nom->GetBinContent(nom->GetMaximumBin());
427  TGraphAsymmErrors* g = new TGraphAsymmErrors;
428  TGraphAsymmErrors* g2 = new TGraphAsymmErrors;
429  for(int binIdx = 0; binIdx < nom->GetNbinsX()+2; ++binIdx){
430  const double y = nom->GetBinContent(binIdx);
431  const double y2 = alt->GetBinContent(binIdx);
432  g->SetPoint(binIdx, nom->GetXaxis()->GetBinCenter(binIdx), y);
433  g2->SetPoint(binIdx, nom->GetXaxis()->GetBinCenter(binIdx), y2);
434  const double w = nom->GetXaxis()->GetBinWidth(binIdx);
435  double errUp = 0, errUp2 = 0;
436  double errDn = 0, errDn2 = 0;
437  for(unsigned int systIdx = 0; systIdx < ups.size(); ++systIdx){
438  double hi = ups[systIdx]->GetBinContent(binIdx)-y;
439  double hi2 = upsA[systIdx]->GetBinContent(binIdx)-y2;
440  double lo = y-dns[systIdx]->GetBinContent(binIdx);
441  double lo2 = y2-dnsA[systIdx]->GetBinContent(binIdx);
442  if(lo <= 0 && hi <= 0) std::swap(lo, hi);
443  if(lo2 <= 0 && hi2 <= 0) std::swap(lo2, hi2);
444  errUp += hi*hi;
445  errDn += lo*lo;
446  errUp2 += hi2*hi2;
447  errDn2 += lo2*lo2;
448  // TODO: what happens if they're both high or both low?
449  } // end for systIdx
450  g->SetPointError(binIdx, w/2, w/2, sqrt(errDn), sqrt(errUp));
451  g2->SetPointError(binIdx, w/2, w/2, sqrt(errDn2), sqrt(errUp2));
452  } // end for i
453  g2->SetFillColorAlpha(col2, 0.3);
454  g2->Draw("2");
455  // g2->GetYaxis()->SetRangeUser(0, headroom*yMax);
456  g2->GetYaxis()->SetRangeUser(0.0, 19.0);
457  g->SetFillColorAlpha(col1, 0.3);
458  g->Draw("2");
459  // g->GetYaxis()->SetRangeUser(0, headroom*yMax);
460  // nom->GetYaxis()->SetRangeUser(0, headroom*yMax);
461  g->GetYaxis()->SetRangeUser(0.0, yMax*headroom);
462  nom->GetYaxis()->SetRangeUser(0.0, yMax*headroom);
463  alt->Draw("hist same");
464  nom->Draw("hist same");
465  for(TH1* up: ups) delete up;
466  for(TH1* dn: dns) delete dn;
467  for(TH1* up2: upsA) delete up2;
468  for(TH1* dn2: dnsA) delete dn2;
469  }
470 
471 
472 
473 
474 
475  //----------------------------------------------------------------------
476  TGraphAsymmErrors* PlotWithSystErrorBand_Quant(const int quant,
477  IPrediction* pred,
478  const std::vector<const ISyst*>& systs,
480  double pot,
481  int col, int errCol,
482  float maxy,
483  bool newaxis)
484  {
485 
486  Spectrum nom = pred->Predict(calc);
487 
488  std::vector<Spectrum> ups, dns;
489 
490  for(const ISyst* syst: systs){
491  SystShifts shifts;
492  shifts.SetShift(syst, +1);
493  ups.push_back(pred->PredictSyst(calc, shifts));
494  shifts.SetShift(syst, -1);
495  dns.push_back(pred->PredictSyst(calc, shifts));
496  }
497 
498  return PlotWithSystErrorBand_Quant(quant, nom, ups, dns, pot, col, errCol, maxy, newaxis);
499 
500 
501  }
502 
503  //----------------------------------------------------------------------
504  TGraphAsymmErrors* PlotWithSystErrorBand_Quant(const int quant,
505  const Spectrum& nominal,
506  const std::vector<Spectrum>& upShifts,
507  const std::vector<Spectrum>& downShifts,
508  double pot, int col, int errCol,
509  float maxy, bool newaxis)
510  {
511 
512  TH1* nom = nominal.ToTH1(pot);
513  nom->GetYaxis()->SetTitle("Events/0.1 GeV");
514  std::vector<TH1*> ups, dns;
515 
516  for(const auto& upShift:upShifts) ups.push_back(upShift.ToTH1(pot));
517  for(const auto& downShift:downShifts) dns.push_back(downShift.ToTH1(pot));
518 
519  return PlotWithSystErrorBand_Quant(quant, nom, ups, dns, col, errCol, maxy, newaxis);
520  }
521 
522  //----------------------------------------------------------------------
523  TGraphAsymmErrors* PlotWithSystErrorBand_Quant(const int quant,
524  TH1*& nom,
525  std::vector<TH1*>& ups,
526  std::vector<TH1*>& dns,
527  int col, int errCol,
528  float maxy, bool newaxis)
529  {
530  if(col == -1){
531  col = kTotalMCColor;
532  errCol = kTotalMCErrorBandColor;
533  }
534  else if(errCol == -1) errCol = col-7; // hopefully a lighter version
535 
536  nom->GetYaxis()->SetTitle("Events/0.1 GeV");
537  nom->SetLineColor(col);
538  nom->GetXaxis()->CenterTitle();
539  nom->GetYaxis()->CenterTitle();
540  if(newaxis) nom->Draw("hist ]["); // Set the axes up
541 
542  TGraphAsymmErrors* g = new TGraphAsymmErrors;
543 
544  for(int binIdx = 0; binIdx < nom->GetNbinsX()+2; ++binIdx){
545  const double y = nom->GetBinContent(binIdx);
546  g->SetPoint(binIdx, nom->GetXaxis()->GetBinCenter(binIdx), y);
547 
548  const double w = nom->GetXaxis()->GetBinWidth(binIdx);
549 
550  double errUp = 0;
551  double errDn = 0;
552 
553  for(unsigned int systIdx = 0; systIdx < ups.size(); ++systIdx){
554  double hi = ups[systIdx]->GetBinContent(binIdx)-y;
555  double lo = dns[systIdx]->GetBinContent(binIdx)-y;
556 
557  // If both shifts are in the same direction use the larger
558  double min = std::min(hi,lo);
559  double max = std::max(hi,lo);
560  if(max < 0) max=0;
561  if(min > 0) min=0;
562 
563  errUp += max*max;
564  errDn += min*min;
565  } // end for systIdx
566 
567  g->SetPointError(binIdx, w/2, w/2, sqrt(errDn), sqrt(errUp));
568  } // end for i
569 
570 
571  // setting 2x2 = 4 pads canvas
572  if(quant==1 || quant==2){
573  nom->GetXaxis()->SetLabelSize(0.);
574  }
575  if(quant==2 || quant==4){
576  nom->GetYaxis()->SetLabelSize(0.);
577  }
578  float minx = 0.99 * nom->GetXaxis()->GetXmin();
579  float maxx = 0.99 * nom->GetXaxis()->GetXmax();
580  nom->GetXaxis()->SetLabelOffset(0.01);
581  nom->GetYaxis()->SetLabelOffset(0.01);
582  nom->GetYaxis()->SetTitleOffset(0.50);
583  nom->GetXaxis()->SetTitleSize(0.0);
584  nom->GetYaxis()->SetTitleSize(0.0);
585  nom->GetXaxis()->SetTitleSize(0.0);
586  nom->GetYaxis()->SetTitleSize(0.0);
587  nom->GetYaxis()->SetRangeUser(0.00001, maxy);
588  nom->GetXaxis()->SetRangeUser(minx, maxx);
589  g->GetYaxis() ->SetRangeUser(0.00001, maxy);
590  g->GetXaxis() ->SetRangeUser(minx, maxx);
591  g->SetFillColor(errCol);
592  g->Draw("e2 same");
593  nom->Draw("hist ][ same");
594 
595  for(TH1* up: ups) delete up;
596  for(TH1* dn: dns) delete dn;
597  return g;
598  }
599 
600 
601 
602 
603 
604 
605  //----------------------------------------------------------------------
606  THStack* ToTHStack(const std::vector<std::pair<const Spectrum&, Color_t>>& s,
607  double pot)
608  {
609  THStack* ret = new THStack;
610  for(auto it: s){
611  TH1* h = it.first.ToTH1(pot);
612  h->SetFillColor(it.second);
613  ret->Add(h, "hist");
614  }
615  return ret;
616  }
617 
618  //----------------------------------------------------------------------
619  /// Helper for \ref AutoPlaceLegend
620  double PointDistanceToBox(double x, double y,
621  double x0, double y0, double x1, double y1)
622  {
623  // Inside
624  if(x > x0 && x < x1 && y > y0 && y < y1) return 0;
625 
626  // Corners
627  double d = util::sqr(x-x0)+util::sqr(y-y0);
628  d = std::min(d, util::sqr(x-x1)+util::sqr(y-y0));
629  d = std::min(d, util::sqr(x-x1)+util::sqr(y-y1));
630  d = std::min(d, util::sqr(x-x0)+util::sqr(y-y1));
631 
632  // Top and bottom edges
633  if(x > x0 && x < x1){
634  d = std::min(d, util::sqr(y-y0));
635  d = std::min(d, util::sqr(y-y1));
636  }
637  // Left and right
638  if(y > y0 && y < y1){
639  d = std::min(d, util::sqr(x-x0));
640  d = std::min(d, util::sqr(x-x1));
641  }
642 
643  return d;
644  }
645 
646  //----------------------------------------------------------------------
647  TLegend* AutoPlaceLegend(double dx, double dy, double yPin)
648  {
649  gPad->Update();
650 
651  // Convert requested width and height into physics coordinates
652  dx *= (gPad->GetX2()-gPad->GetX1());
653  dy *= (gPad->GetY2()-gPad->GetY1());
654 
655  // Range of axes in physics units
656  const double x0 = gPad->GetUxmin();
657  const double x1 = gPad->GetUxmax();
658  const double y0 = gPad->GetUymin();
659  const double y1 = gPad->GetUymax();
660 
661  const double X = x1-x0;
662  const double Y = y1-y0;
663 
664  double bestd = 0;
665  double bestx = 0;
666  double besty = 0;
667  // If we want to pin Y pos, set it to that now.
668  if(yPin >= 0) besty = yPin * (gPad->GetY2()-gPad->GetY1());
669 
670  for(bool fallback: {false, true}){
671  for(double x = x0+dx/2; x < x1-dx/2; x += X/50){
672  for(double y = y0+dy/2; y < y1-dy/2; y += Y/50){
673  double d = 999999;
674 
675  // Repel from edges
676  d = std::min(d, util::sqr((x-dx/2-x0)/X));
677  d = std::min(d, util::sqr((x+dx/2-x1)/X));
678  d = std::min(d, util::sqr((y-dy/2-y0)/Y));
679  d = std::min(d, util::sqr((y+dy/2-y1)/Y));
680  if(d < bestd) continue;
681 
682  TIter next(gPad->GetListOfPrimitives());
683  while(TObject* obj = next()){
684 
685  if(!obj->InheritsFrom(TH1::Class())) continue;
686  if( obj->InheritsFrom(TH2::Class())) continue;
687 
688  TH1* h = (TH1*)obj;
689 
690  for(int n = 1; n <= h->GetNbinsX(); ++n){
691  const double px = h->GetBinCenter(n);
692  const double py = h->GetBinContent(n);
693 
694  if(fallback){
695  d = std::min(d, util::sqr((px-x)/X)+util::sqr((py-y)/Y));
696  }
697  else{
698  d = std::min(d, PointDistanceToBox(px/X, py/Y,
699  (x-dx/2)/X, (y-dy/2)/Y,
700  (x+dx/2)/X, (y+dy/2)/Y));
701  }
702  if(d < bestd) break;
703  }
704  if(d < bestd) break;
705  } // end while
706 
707  if(d > bestd){
708  bestd = d;
709  bestx = x;
710  // Update Y if we're not pinning it.
711  if (yPin < 0) besty = y;
712  }
713  } // end for y
714  } // end for x
715 
716  if(bestd != 0) break; // If we always collide, have to do fallback
717  } // end for fallback
718 
719  // Convert to pad coordinates
720  const double nx = (bestx-gPad->GetX1())/(gPad->GetX2()-gPad->GetX1());
721  const double ny = (besty-gPad->GetY1())/(gPad->GetY2()-gPad->GetY1());
722 
723  const double ndx = dx/(gPad->GetX2()-gPad->GetX1());
724  const double ndy = dy/(gPad->GetY2()-gPad->GetY1());
725 
726  return new TLegend(nx-ndx/2, ny-ndy/2, nx+ndx/2, ny+ndy/2);
727  }
728 
729  //----------------------------------------------------------------------
730  void ErrorBarChart(const std::map<std::string, std::pair<double, double>>& systErrors,
731  const std::pair<double, double> & statErr,
732  const std::string & label )
733  {
734  bool doStat = statErr.first != 0 && statErr.second != 0;
735  int nbins = systErrors.size() + 1; // extra for the "total systematic" bar
736 
737  // Systematics are on top, total systematic and (optional) statistical
738  // error bar are at the bottom.
739  int firstSystBin = 1;
740  if(doStat){
741  ++nbins;
742  ++firstSystBin;
743  }
744 
745  // Sum all the systematics in qaudrature
746  std::pair<double, double> systTot = std::make_pair(0., 0.);
747  for(auto it_systPairs: systErrors)
748  {
749  systTot.first += util::sqr(it_systPairs.second.first);
750  systTot.second += util::sqr(it_systPairs.second.second);
751  }
752  systTot.first = sqrt(systTot.first);
753  systTot.second = sqrt(systTot.second);
754 
755  // x-axis range is a multiple of 5% with at least 2% padding
756  double xrange = 0;
757  while(xrange < std::max(std::max(systTot.first, systTot.second), std::max(statErr.first, statErr.second))+2) xrange += 5;
758 
759  TH2* axes = new TH2F("", (";" + label).c_str(),
760  10, -xrange, +xrange, nbins, 0, nbins);
761  axes->GetXaxis()->CenterTitle();
762 
763  // Put the systematics in a vector so we can sort them by size
764  std::vector<std::pair<std::pair<double, double>, std::string>> systList;
765  for(auto it: systErrors) systList.push_back(std::make_pair(it.second, it.first));
766  std::sort(systList.begin(),
767  systList.end(),
768  [](const auto & errNamePair1, const auto & errNamePair2)
769  {
770  return ( std::max(errNamePair1.first.first, errNamePair1.first.second)
771  < std::max(errNamePair2.first.first, errNamePair2.first.second) );
772  }
773  );
774 
775  // Label the y axis
776  TAxis* yax = axes->GetYaxis();
777  for(unsigned int i = 0; i < systList.size(); ++i)
778  yax->SetBinLabel(i+firstSystBin+1, systList[i].second.c_str());
779 
780  yax->SetBinLabel(firstSystBin, "Total syst. error");
781  if(doStat) yax->SetBinLabel(1, "Statistical error");
782 
783  // Make the labels approximately match the x axis title
784  yax->SetLabelSize(.06);
785 
786  axes->Draw();
787 
788  // Blue boxes for each systematic
789  for(unsigned int i = 0; i < systList.size(); ++i){
790  TBox* box = new TBox(-systList[i].first.first, i+firstSystBin+.2,
791  +systList[i].first.second, i+firstSystBin+.8);
792  std::cout<<"ERROR = (" << systList[i].first.first << ", " << systList[i].first.second << ")" << std::endl;
793  box->SetFillColor(kAzure-8);
794  box->Draw();
795  }
796 
797  // Total systematic error
798  TBox* syst = new TBox(-systTot.first, firstSystBin-.8,
799  +systTot.second, firstSystBin-.2);
800  std::cout<<"TOTAL SYSTEMATIC ERROR = (" << systTot.first << ", " << systTot.second << ")" << std::endl;
801  syst->SetFillColor(kAzure-8);
802  syst->Draw();
803 
804  // Red box for statistical error
805  if(doStat){
806  TBox* stat = new TBox(-statErr.first, .2, +statErr.second, .8);
807  stat->SetFillColor(kRed-7);
808  stat->Draw();
809  }
810 
811  // Separate individual systematics from the total and statistical errors
812  TLine* div = new TLine(-xrange, firstSystBin, +xrange, firstSystBin);
813  div->SetLineWidth(2);
814  div->Draw();
815 
816  // Vertical line marking the symmetry point
817  TLine* zero = new TLine(0, 0, 0, nbins);
818  zero->SetLineStyle(7);
819  zero->SetLineWidth(2);
820  zero->Draw();
821 
822  // Leave enough space for the systematic labels
823  gPad->SetLeftMargin(.25);
824 
825  }
826 
827  //----------------------------------------------------------------------
828  void CountingExperimentErrorBarChart(const std::map<std::string, double>& systbars, double statErr, bool bkgdOrSig, bool shortchart)
829  {
830  std::map<std::string, double> systs = (shortchart) ? SumNueSecondAnaSysts(systbars) : systbars;
831 
832  std::map<std::string, std::pair<double, double>> systVals;
833  for (const auto & systPair : systbars)
834  systVals[systPair.first] = std::make_pair(systPair.second, systPair.second);
835  auto statVals = std::make_pair(statErr, statErr);
836 
837  std::string label = std::string(bkgdOrSig ? "Signal" : "Background") + " uncertainty (%)";
838  ErrorBarChart(systVals, statVals, label);
839  }
840 
841  //----------------------------------------------------------------------
842  TGraphAsymmErrors* GraphWithPoissonErrors(const TH1* h, bool noErrorsXaxis, bool drawEmptyBins)
843  {
844  TGraphAsymmErrors* gr = new TGraphAsymmErrors(h);
845 
846  TFeldmanCousins fc(0.6827);//1 sigma
847 
848  for(int i = 0; i < h->GetNbinsX(); ++i){
849  double x, y;
850  gr->GetPoint(i, x, y);
851 
852  if ( drawEmptyBins || y!=0 ){
853  if ( y < 50 ) gr->SetPointEYlow(i, y-fc.CalculateLowerLimit(y,0));
854  else gr->SetPointEYlow(i,sqrt(y));
855  if ( y < 30 ) gr->SetPointEYhigh(i,fc.CalculateUpperLimit(y,0)-y);
856  else gr->SetPointEYhigh(i,sqrt(y));
857  }
858 
859  if(noErrorsXaxis){
860  gr->SetPointEXlow(i, 0);
861  gr->SetPointEXhigh(i, 0);
862  } // Do not use bin width as X-error for points
863  }
864 
865  return gr;
866  }
867 
868 
869  //----------------------------------------------------------------------
870  TGraph* graphAsymmErrorScaled(TH1* histScaled, TH1* hist, double overallScale){
871 
872  for(int bin = 1; bin <= hist->GetNbinsX(); bin++) {
873  if(hist->GetBinContent(bin)==0){
874  hist->SetBinContent(bin,0.001);
875  histScaled->SetBinContent(bin,0.001);
876  }
877  }
878 
879  TFeldmanCousins fc(0.6827); // to get correct poisson error bars
880  TGraphAsymmErrors *datapoisson = new TGraphAsymmErrors(histScaled);
881  for(int bin = 1; bin <= hist->GetNbinsX(); bin++) {
882  //default numu: scale by bin 0.1/width
883  float binWidth = hist->GetBinWidth(bin);
884  float scaleFactor = overallScale / binWidth;
885  float y = hist->GetBinContent(bin);
886  float errup = 0, errdn = 0;
887  if(y < 50) errdn = scaleFactor * (y - fc.CalculateLowerLimit(y,0));
888  else errdn = scaleFactor * sqrt(y);
889  if(y < 30) errup = scaleFactor * (fc.CalculateUpperLimit(y,0) - y);
890  else errup = scaleFactor * sqrt(y);
891 
892  datapoisson->SetPointEYhigh(bin-1, errup);
893  datapoisson->SetPointEYlow(bin-1, errdn);
894  }
895  datapoisson->SetMarkerStyle(kFullCircle);
896  datapoisson->SetMarkerColor(1);
897  datapoisson->SetMarkerSize(1);
898  datapoisson->SetLineWidth(2);
899 
900  return datapoisson;
901 
902  }
903 
904 
905 
906  //----------------------------------------------------------------------
907  TGraph* ShadeBetweenHistograms(TH1* hmin, TH1* hmax)
908  {
909  int n = hmin->GetNbinsX();
910  TGraph* gr = new TGraph(4*n);
911 
912  for(int i=1; i<=n; i++)
913  {
914  double xdown = hmax->GetBinLowEdge(i);
915  double xup = hmax->GetBinLowEdge(i+1);
916  double ydown = hmin->GetBinContent(i);
917  double yup = hmax->GetBinContent(i);
918 
919  gr->SetPoint(2*(i-1), xdown, ydown);
920  gr->SetPoint(2*(i-1)+1, xup, ydown);
921  gr->SetPoint(4*n-2*i, xup, yup);
922  gr->SetPoint(4*n-2*i+1, xdown, yup);
923  }
924 
925  gr->SetFillColor(hmin->GetLineColor() - 9); // In principle this is lighter
926  gr->SetLineColor(hmin->GetLineColor() - 9); // In case one does Draw("lf")
927 
928  return gr;
929  }
930  //----------------------------------------------------------------------
931  TGraphAsymmErrors * ProfileQuantile(const TH2 * hist,
932  const std::string & axisName,
933  const std::string & graphName,
934  const std::pair<double, double> & quantileDivisions)
935  {
936  if (hist->GetDimension() != 2)
937  throw std::runtime_error(Form("Can't profile a histogram with other than 2 dimensions. Yours is a %d-D histogram...", hist->GetDimension()));
938 
939  TAxis const* axis = nullptr;
940  std::function<TH1D*(const char *, Int_t, Int_t, Option_t*)> projectionMethod;
941  if (axisName == "x" || axisName == "X")
942  {
943  axis = hist->GetXaxis();
944  projectionMethod = std::bind(&TH2::ProjectionY, hist, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4);
945  }
946  else if (axisName == "y" || axisName == "Y")
947  {
948  axis = hist->GetYaxis();
949  projectionMethod = std::bind(&TH2::ProjectionX, hist, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4);
950  }
951  else
952  throw std::runtime_error( Form("Can't profile onto unknown axis: '%s'", axisName.c_str()) );
953 
954 
955  std::vector<double> points_x, points_y, errors_x_low, errors_x_high, errors_y_low, errors_y_high;
956 
957  double quantiles[2] = {0, 0};
958  double _quantileDivisions[2] = {quantileDivisions.first, quantileDivisions.second};
959  for (int binNum = 1; binNum <= axis->GetNbins(); binNum++) // remember, 0 is underflow and Nbins+1 is overflow...
960  {
961  std::unique_ptr<TH1D> proj ( projectionMethod(Form("%s_tmp_%d", hist->GetName(), int(std::time(nullptr))), binNum, binNum, "") ); // randomish name, just to be safe
962 
963  double y = proj->GetMean();
964  points_x.push_back(axis->GetBinCenter(binNum));
965  points_y.push_back(y);
966 
967  // for now, make the errors for an empty projection 0
968  unsigned int n_qs = 0;
969  if (proj->Integral() == 0)
970  {
971  n_qs = 2;
972  quantiles[0] = quantiles[1] = 0;
973  }
974  else
975  n_qs = proj->GetQuantiles(2, quantiles, _quantileDivisions);
976  if (n_qs != 2)
977  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) );
978 
979  double binWidth = axis->GetBinWidth(binNum);
980  errors_x_low.push_back(binWidth/2.);
981  errors_x_high.push_back(binWidth/2.);
982  errors_y_low.push_back(y-quantiles[0]);
983  errors_y_high.push_back(quantiles[1]-y);
984  }
985 
986  TGraphAsymmErrors * outGraph = new TGraphAsymmErrors(
987  points_x.size(),
988  &points_x[0],
989  &points_y[0],
990  &errors_x_low[0],
991  &errors_x_high[0],
992  &errors_y_low[0],
993  &errors_y_high[0]
994  );
995  std::string name = (graphName.size()) ? graphName : Form("%s_quantile_errors", hist->GetName());
996  outGraph->SetName(name.c_str());
997 
998  return outGraph;
999 
1000  }
1001 
1002 
1003  // Draw single BF
1004  void drawBFSingle(double bfSin, double bfDm, Color_t color, Style_t marker, double size = 1.5)
1005  {
1006 
1007  TMarker* manMarker = new TMarker(bfSin, bfDm, marker);
1008  manMarker->SetMarkerSize(size);
1009  manMarker->SetMarkerColor(color);
1010  manMarker->Draw();
1011 
1012  }
1013 
1014  // Truth for FHC; not sure about RHC
1015  void drawBFMirror(double bfSin, double bfDm, Color_t color, Style_t marker, double size = 1.5)
1016  {
1017 
1018  double mirror = bfSin;
1019  TMarker* manMarker = new TMarker(bfSin, bfDm, marker);
1020  manMarker->SetMarkerSize(size);
1021  manMarker->SetMarkerColor(color);
1022  manMarker->Draw();
1023 
1024  if(bfSin > 0.514) mirror = 0.514 - (bfSin - 0.514);
1025  if(bfSin < 0.514) mirror = 0.514 + (0.514 - bfSin );
1026  TMarker* mirrorMarker = new TMarker(mirror, bfDm, marker);
1027  mirrorMarker->SetMarkerSize(size);
1028  mirrorMarker->SetMarkerColor(color);
1029  mirrorMarker->Draw();
1030 
1031  }
1032 
1033 
1034 
1035  //----------------------------------------------------------------------
1036  void MakeHistCanvasReady_Quant(const int quant,
1037  TH1* hist,
1038  float maxy){
1039 
1040  if(quant==1 || quant==2) hist->GetXaxis()->SetLabelSize(0.);
1041  if(quant==2 || quant==4) hist->GetYaxis()->SetLabelSize(0.);
1042  if(quant == 3) hist->GetXaxis()->ChangeLabel(6,-1,0);
1043 
1044  hist->GetXaxis()->SetLabelOffset(0.01);
1045  hist->GetYaxis()->SetLabelOffset(0.01);
1046  hist->GetYaxis()->SetTitleOffset(0.50);
1047  hist->GetXaxis()->SetTitleSize(0.0);
1048  hist->GetYaxis()->SetTitleSize(0.0);
1049  hist->GetYaxis()->SetRangeUser(0.0, maxy);
1050  }
1051 
1052  //----------------------------------------------------------------------
1053  void PimpHist(TH1* hist, Style_t linestyle, Color_t linecolor, int linewidth, Style_t markerstyle, Color_t markercolor, double markersize){
1054  hist->SetLineColor(linecolor);
1055  hist->SetLineStyle(linestyle);
1056  hist->SetLineWidth(linewidth);
1057  hist->SetMarkerColor(markercolor);
1058  hist->SetMarkerStyle(markerstyle);
1059  hist->SetMarkerSize(markersize);
1060  }
1061 
1062 
1063  //----------------------------------------------------------------------
1064  void SplitCanvas(double ysplit, TPad*& p1, TPad*& p2)
1065  {
1066  if(!gPad) new TCanvas;
1067  p1 = new TPad("", "", 0, 0, 1, 1);
1068  p2 = new TPad("", "", 0, 0, 1, 1);
1069 
1070  p1->SetBottomMargin(ysplit);
1071  p2->SetTopMargin(1-ysplit);
1072 
1073  // Draw p1 second since it's often the more important one, that the user
1074  // would prefer to be able to interact with.
1075  for(TPad* p: {p2, p1}){
1076  p->SetFillStyle(0);
1077  p->Draw();
1078  }
1079  }
1080 
1081 
1082  //----------------------------------------------------------------------
1083  // split canvas in 4 for numu quartiles
1084  void SplitCanvasQuant(TCanvas *& canvas, TPad *& pad1, TPad *& pad2, TPad *& pad3, TPad *& pad4){
1085 
1086  canvas = new TCanvas(UniqueName().c_str(),"",960,800);
1087 
1088  pad3 = new TPad(UniqueName().c_str(),"pad3",0,0,1,1);
1089  pad3->Draw();
1090  pad3->SetTopMargin(0.5);
1091  pad3->SetBottomMargin(0.1);
1092  pad3->SetLeftMargin(0.08);
1093  pad3->SetRightMargin(0.5);
1094  pad3->SetFillStyle(0);
1095  canvas->cd();
1096 
1097  pad4 = new TPad(UniqueName().c_str(),"pad4",0,0,1,1);
1098  pad4->Draw();
1099  pad4->SetTopMargin(0.5);
1100  pad4->SetBottomMargin(0.1);
1101  pad4->SetLeftMargin(0.5);
1102  pad4->SetRightMargin(0.08);
1103  pad4->SetFillStyle(0);
1104  canvas->cd();
1105 
1106  pad1 = new TPad(UniqueName().c_str(),"pad1",0,0,1,1);// x1 y1 x2 y2
1107  pad1->Draw();
1108  pad1->SetTopMargin(0.1);
1109  pad1->SetBottomMargin(0.5);
1110  pad1->SetLeftMargin(0.08);
1111  pad1->SetRightMargin(0.5);
1112  pad1->SetFillStyle(0);
1113  canvas->cd();
1114 
1115  pad2 = new TPad(UniqueName().c_str(),"pad2",0,0,1,1);
1116  pad2->Draw();
1117  pad2->SetTopMargin(0.1);
1118  pad2->SetBottomMargin(0.5);
1119  pad2->SetLeftMargin(0.5);
1120  pad2->SetRightMargin(0.08);
1121  pad2->SetFillStyle(0);
1122  canvas->cd();
1123 
1124  }
1125 
1126 
1127  //----------------------------------------------------------------------
1128  void CenterTitles ( TH1 * histo)
1129  {
1130  histo->GetXaxis()->CenterTitle();
1131  histo->GetYaxis()->CenterTitle();
1132  histo->GetZaxis()->CenterTitle();
1133  }
1134 
1135 
1136  //----------------------------------------------------------------------
1137  bool SortSystsName(const ISyst* systA, const ISyst* systB){
1138  return ( systA->ShortName() < systB->ShortName() );
1139  }
1140 
1141 
1142  TH1* PlotSystShifts(const SystShifts & shifts, bool sortName){
1143 
1144  auto systs = shifts.ActiveSysts();
1145  int nsysts = systs.size();
1146 
1147  // sort systematics by name
1148  if(sortName) std::sort(systs.begin(), systs.end(), SortSystsName);
1149 
1150  auto h = new TH1D ("sh",";; Pull (N#sigma)",
1151  nsysts, -0.5, nsysts + 0.5);
1152 
1153  for (int systIdx = 0; systIdx < nsysts; ++systIdx){
1154  double shiftThis = shifts.GetShift(systs[systIdx]);
1155  h->SetBinContent( systIdx + 1, shiftThis);
1156 
1157  TString tempName = systs[systIdx]->ShortName();
1158  if( tempName.Contains("Neutron") ) h->GetXaxis()->SetBinLabel( systIdx + 1, "NeutronSyst2018");
1159  else h->GetXaxis()->SetBinLabel( systIdx + 1, systs[systIdx]->ShortName().c_str());
1160  }
1161 
1162  h->SetMarkerStyle(kFullCircle);
1163  CenterTitles(h);
1164  h->SetTitleOffset(3.0);
1165  h->GetXaxis()->SetNdivisions(nsysts,kFALSE);
1166  h->GetXaxis()->SetLabelSize(0.065);
1167  h->GetXaxis()->LabelsOption("v");
1168  h->GetYaxis()->SetRangeUser(-1,1);
1169  auto c1 = new TCanvas ("c1","c1",700,350);
1170  c1->SetBottomMargin(0.5);
1171  h->Draw("lp hist");
1172  return h;
1173  }
1174 
1175  //----------------------------------------------------------------------
1176  TGraph* JoinGraphs(TGraph* a, TGraph* b, int fillcol)
1177  {
1178  TGraph* ret = (TGraph*)a->Clone(UniqueName().c_str());
1179  for(int i = 0; i < b->GetN(); ++i){
1180  double x, y;
1181  b->GetPoint(i, x, y);
1182  ret->SetPoint(ret->GetN(), x, y);
1183  }
1184 
1185  ret->SetFillColor(fillcol);
1186 
1187  return ret;
1188  }
1189 
1190  //----------------------------------------------------------------------
1191 
1192  TGraph* ExtendGraphToTop(TGraph* g, int col, double xmin, double xmax, double y)
1193  {
1194  g->SetPoint(g->GetN(), xmax, y);
1195  g->SetPoint(g->GetN(), xmin, y);
1196 
1197  g->SetFillColor(col);
1198 
1199  return g;
1200  }
1201 
1202 
1203 } // namespace
T max(const caf::Proxy< T > &a, T b)
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")
osc::OscCalculatorDumb calc
Oscillation analysis framework, runs over CAF files outside of ART.
set< int >::iterator it
TH1 * PlotSystShifts(const SystShifts &shifts, bool sortName)
Definition: Plots.cxx:1142
std::vector< SystGroupDef > systs
Definition: syst_header.h:384
TH1D * ToTH1(Color_t col=kBlack, Style_t style=kSolid) const
Definition: Ratio.cxx:122
General interface to oscillation calculators.
Definition: FwdDeclare.h:15
void MakeHistCanvasReady_Quant(const int quant, TH1 *hist, float maxy)
Definition: Plots.cxx:1036
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:553
virtual Spectrum PredictComponent(osc::IOscCalculator *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
virtual const std::string & ShortName() const final
The name printed out to the screen.
Definition: ISyst.h:27
TH1 * DataMCComparisonComponents(const Spectrum &data, const IPrediction *mc, osc::IOscCalculator *calc)
Plot MC broken down into flavour components, overlayed with data.
Definition: Plots.cxx:113
void DataMCRatio(const Spectrum &data, const IPrediction *mc, osc::IOscCalculator *calc, double miny, double maxy)
Plot data/MC ratio for the given spectrum. Normalize MC to Data by POT.
Definition: Plots.cxx:166
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:269
double maxy
Float_t x1[n_points_granero]
Definition: compare.C:5
(&#39;beam &#39;)
Definition: IPrediction.h:15
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
TGraphAsymmErrors * PlotWithSystErrorBand(IPrediction *pred, const std::vector< const ISyst * > &systs, osc::IOscCalculator *calc, double pot, int col, int errCol, float headroom, bool newaxis, EBinType bintype)
Plot prediction with +/-1sigma error band.
Definition: Plots.cxx:244
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:335
void CenterTitles(TH1 *histo)
Definition: Plots.cxx:1128
TGraphAsymmErrors * GraphWithPoissonErrors(const TH1 *h, bool noErrorsXaxis, bool drawEmptyBins)
Calculate statistical errors appropriate for small Poisson numbers.
Definition: Plots.cxx:842
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:1064
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:742
const Color_t kTotalMCErrorBandColor
Definition: Style.h:18
::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
std::map< std::string, double > SumNueSecondAnaSysts(std::map< std::string, double > systbars)
const Color_t kNumuBackgroundColor
Definition: Style.h:31
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:931
virtual Spectrum Predict(osc::IOscCalculator *calc) const =0
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:37
virtual Spectrum PredictSyst(osc::IOscCalculator *calc, const SystShifts &syst) const
Definition: IPrediction.cxx:98
Float_t Y
Definition: plot.C:38
void PlotWithAreaSystErrorBand(IPrediction *pred, const std::vector< const ISyst * > &systs, osc::IOscCalculator *calc, double pot, int col, int errCol, float headroom, bool newaxis, EBinType bintype)
Plot prediction with +/-1sigma shape-only error band.
Definition: Plots.cxx:352
const char * label
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:1053
const XML_Char const XML_Char * data
Definition: expat.h:268
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
TGraphAsymmErrors * PlotWithSystErrorBand_Quant(const int quant, IPrediction *pred, const std::vector< const ISyst * > &systs, osc::IOscCalculator *calc, double pot, int col, int errCol, float maxy, bool newaxis)
Definition: Plots.cxx:476
EBinType
Definition: Utilities.h:46
void drawBFSingle(double bfSin, double bfDm, Color_t color, Style_t marker, double size=1.5)
Draw a single BF point.
Definition: Plots.cxx:1004
const XML_Char * s
Definition: expat.h:262
const int nbins
Definition: cellShifts.C:15
Charged-current interactions.
Definition: IPrediction.h:39
TSpline3 hi("hi", xhi, yhi, 18,"0")
double dy[NP][NC]
double dx[NP][NC]
void RatioPlot(const Spectrum &data, const Spectrum &expected, const Spectrum &fit, double miny, double maxy)
Plot data/expected, compared with fit/expected.
Definition: Plots.cxx:215
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:403
TH1 * DataMCComparisonAreaNormalized(const Spectrum &data, const Spectrum &mc)
Definition: Plots.cxx:73
TGraph * graphAsymmErrorScaled(TH1 *histScaled, TH1 *hist, double overallScale)
Definition: Plots.cxx:870
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:907
#define pot
T GetShift(const ISyst *syst) const
Float_t d
Definition: plot.C:236
TGraph * ExtendGraphToTop(TGraph *g, int col, double xmin, double xmax, double y)
Definition: Plots.cxx:1192
TH2D * histo
float bin[41]
Definition: plottest35.C:14
double POT() const
Definition: Spectrum.h:263
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:1176
OStream cout
Definition: OStream.cxx:6
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:620
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:828
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void DataMCAreaNormalizedRatio(const Spectrum &data, const IPrediction *mc, osc::IOscCalculator *calc, double miny, double maxy)
Plot data/MC ratio for the given spectrum. Normalize MC to Data by area.
Definition: Plots.cxx:195
string syst
Definition: plotSysts.py:176
Float_t proj
Definition: plot.C:35
const hit & b
Definition: hits.cxx:21
void ProjectionY(T *from, U *to)
Helper for WeightingVariable.
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:730
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:606
bool SortSystsName(const ISyst *systA, const ISyst *systB)
Definition: Plots.cxx:1137
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
c1
Definition: demo5.py:24
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
TPad * pad3
Definition: analysis.C:13
std::vector< const ISyst * > ActiveSysts() const
Definition: SystShifts.cxx:199
string tempName
TLegend * AutoPlaceLegend(double dx, double dy, double yPin)
Create a legend, maximizing distance from all histograms.
Definition: Plots.cxx:647
const Color_t kTotalMCColor
Definition: Style.h:17
void SplitCanvasQuant(TCanvas *&canvas, TPad *&pad1, TPad *&pad2, TPad *&pad3, TPad *&pad4)
Definition: Plots.cxx:1084
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:34
T min(const caf::Proxy< T > &a, T b)
All neutrinos, any flavor.
Definition: IPrediction.h:26
Float_t X
Definition: plot.C:38
const Color_t kNCBackgroundColor
Definition: Style.h:22
Float_t w
Definition: plot.C:20
void ProjectionX(T *from, U *to)
Helper for Unweighted.
void next()
Definition: show_event.C:84
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:70
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:1015
TPad * pad1
Definition: analysis.C:13
h
Definition: demo3.py:41