All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Plots.cxx
Go to the documentation of this file.
1 #include "sbnana/CAFAna/Analysis/Plots.h"
2 
3 #include "sbnana/CAFAna/Analysis/Style.h"
4 #include "sbnana/CAFAna/Prediction/IPrediction.h"
5 #include "sbnana/CAFAna/Core/Ratio.h"
6 #include "sbnana/CAFAna/Core/Spectrum.h"
7 #include "sbnana/CAFAna/Core/SystShifts.h"
8 
9 #include "sbnana/CAFAna/Core/MathUtil.h"
10 
11 #include "TCanvas.h"
12 #include "TFeldmanCousins.h"
13 #include "TGraph.h"
14 #include "TGraphAsymmErrors.h"
15 #include "TH1.h"
16 #include "TH2.h"
17 #include "THStack.h"
18 #include "TLegend.h"
19 #include "TLine.h"
20 #include "TList.h"
21 
22 #include <algorithm>
23 #include <iostream>
24 #include <cmath>
25 #include <ctime>
26 #include <functional>
27 #include <memory>
28 
29 namespace ana
30 {
31  //----------------------------------------------------------------------
32  TH1* DataMCComparison(const Spectrum& data, const Spectrum& mc)
33  {
34  TH1* ret = 0;
35 
36  TH1* hMC = mc.ToTH1(data.POT());
37  hMC->SetLineColor(kTotalMCColor);
38 
39  TH1* hData = data.ToTH1(data.POT());
40  hData->Sumw2();
41  hData->SetMarkerStyle(kFullCircle);
42 
43  // Need to draw the biggest thing first to get correct axis limits
44  if(hMC->GetMaximum() > hData->GetMaximum()+sqrt(hData->GetMaximum())){
45  ret = hMC;
46  hMC->Draw("hist");
47  hData->Draw("ep same");
48  }
49  else{
50  ret = hData;
51  hData->Draw("ep");
52  hMC->Draw("hist same");
53  hData->Draw("ep same"); // Data always goes on top
54  }
55 
56  gPad->Update();
57 
58  return ret;
59  }
60 
61  //----------------------------------------------------------------------
62  TH1* DataMCComparisonAreaNormalized(const Spectrum& data, const Spectrum& mc)
63  {
64  TH1* ret = 0;
65 
66  TH1* hMC = mc.ToTH1(mc.POT());
67  hMC->SetLineColor(kTotalMCColor);
68 
69  TH1* hData = data.ToTH1(data.POT());
70  hData->Sumw2();
71  hData->SetMarkerStyle(kFullCircle);
72 
73  hMC->Scale(hData->Integral()/hMC->Integral());
74 
75  // Need to draw the biggest thing first to get correct axis limits
76  if(hMC->GetMaximum() > hData->GetMaximum()+sqrt(hData->GetMaximum())){
77  ret = hMC;
78  hMC->Draw("hist");
79  hData->Draw("ep same");
80  }
81  else{
82  ret = hData;
83  hData->Draw("ep");
84  hMC->Draw("hist same");
85  hData->Draw("ep same"); // Data always goes on top
86  }
87 
88  gPad->Update();
89 
90  return ret;
91  }
92 
93  //----------------------------------------------------------------------
94  TH1* DataMCComparisonComponents(const Spectrum& data,
95  const IPrediction* mc,
96  osc::IOscCalc* calc)
97  {
98  TH1* ret = 0;
99 
100  TH1* hMC = mc->Predict(calc).ToTH1(data.POT());
101  hMC->SetLineColor(kTotalMCColor);
102 
103  TH1* hNC = mc->PredictComponent(calc,
104  Flavors::kAll,
105  Current::kNC,
106  Sign::kBoth).ToTH1(data.POT());
107  hNC->SetLineColor(kNCBackgroundColor);
108 
109  TH1* hCC = mc->PredictComponent(calc,
110  Flavors::kAllNuMu,
111  Current::kCC,
112  Sign::kBoth).ToTH1(data.POT());
113  hCC->SetLineColor(kNumuBackgroundColor);
114 
115  TH1* hBeam = mc->PredictComponent(calc,
116  Flavors::kNuEToNuE,
117  Current::kCC,
118  Sign::kBoth).ToTH1(data.POT());
119  hBeam->SetLineColor(kBeamNueBackgroundColor);
120 
121  TH1* hData = data.ToTH1(data.POT());
122  hData->SetMarkerStyle(kFullCircle);
123 
124  // Need to draw the biggest thing first to get correct axis limits
125  if(hMC->GetMaximum() > hData->GetMaximum()+sqrt(hData->GetMaximum())){
126  ret = hMC;
127  hMC->Draw("hist");
128  hData->Draw("ep same");
129  }
130  else{
131  ret = hData;
132  hData->Draw("ep");
133  hMC->Draw("hist same");
134  hData->Draw("ep same"); // Data always goes on top
135  }
136 
137  hNC->Draw("hist same");
138  hCC->Draw("hist same");
139  hBeam->Draw("hist same");
140 
141  gPad->Update();
142 
143  return ret;
144  }
145 
146  TH1* GetMCSystTotal(const IPrediction* mc,
147  osc::IOscCalc* calc,
148  const SystShifts& shift,
149  std::string hist_name,
150  double pot,
151  bool force1D)
152  {
153  TH1* hTotal = mc->PredictSyst(calc, shift).ToTHX(pot, force1D);
154  hTotal->SetNameTitle((hist_name+"_total").c_str(), (hist_name+"_total").c_str());
155  return hTotal;
156  }
157 
158  TH1* GetMCTotal(const IPrediction* mc,
159  osc::IOscCalc* calc,
160  std::string hist_name,
161  double pot,
162  bool force1D)
163  {
164  return GetMCSystTotal(mc, calc, kNoShift, hist_name, pot, force1D);
165  }
166 
167  std::vector<TH1*> GetMCComponents(const IPrediction* mc,
168  osc::IOscCalc* calc,
169  std::string hist_name,
170  double pot,
171  bool force1D)
172  {
173  return GetMCSystComponents(mc, calc, kNoShift, hist_name, pot, force1D);
174  }
175 
176  std::vector<TH1*> GetMCSystComponents(const IPrediction* mc,
177  osc::IOscCalc* calc,
178  const SystShifts& shift,
179  std::string hist_name,
180  double pot,
181  bool force1D)
182  {
183  std::vector<TH1*> ret;
184 
185  TH1* hTotal = mc->PredictSyst(calc, shift).ToTHX(pot, force1D);
186  hTotal->SetNameTitle((hist_name+"_total").c_str(), (hist_name+"_total").c_str());
187  ret .push_back(hTotal);
188 
189  TH1* hNC = mc->PredictComponentSyst(calc,shift,
190  Flavors::kAll,
191  Current::kNC,
192  Sign::kBoth).ToTHX(pot, force1D);
193  hNC->SetNameTitle((hist_name+"_NC").c_str(), (hist_name+"_NC").c_str());
194  ret .push_back(hNC);
195 
196  TH1* hAllNumu = mc->PredictComponentSyst(calc,shift,
197  Flavors::kAllNuMu,
198  Current::kCC,
199  Sign::kBoth).ToTHX(pot, force1D);
200  hAllNumu->SetNameTitle((hist_name+"_AllNumu").c_str(), (hist_name+"_AllNumu").c_str());
201  ret .push_back(hAllNumu);
202 
203  TH1* hNumu = mc->PredictComponentSyst(calc,shift,
204  Flavors::kAllNuMu,
205  Current::kCC,
206  Sign::kNu).ToTHX(pot, force1D);
207  hNumu->SetNameTitle((hist_name+"_Numu").c_str(), (hist_name+"_Numu").c_str());
208  ret .push_back(hNumu);
209 
210  TH1* hNumubar = mc->PredictComponentSyst(calc,shift,
211  Flavors::kAllNuMu,
212  Current::kCC,
213  Sign::kAntiNu).ToTHX(pot, force1D);
214  hNumubar->SetNameTitle((hist_name+"_Numubar").c_str(), (hist_name+"_Numubar").c_str());
215  ret .push_back(hNumubar);
216 
217  TH1* hAllNutau = mc->PredictComponentSyst(calc,shift,
218  Flavors::kAllNuTau,
219  Current::kCC,
220  Sign::kBoth).ToTHX(pot, force1D);
221  hAllNutau->SetNameTitle((hist_name+"_AllNutau").c_str(), (hist_name+"_AllNutau").c_str());
222  ret .push_back(hAllNutau);
223 
224  // Want AllSignalNue, SignalNue, SignalNuebar
225  TH1* hAllNue = mc->PredictComponentSyst(calc,shift,
226  Flavors::kAllNuE,
227  Current::kCC,
228  Sign::kBoth).ToTHX(pot, force1D);
229  hAllNue->SetNameTitle((hist_name+"_AllNue").c_str(), (hist_name+"_AllNue").c_str());
230  ret .push_back(hAllNue);
231 
232  TH1* hAllBeamNue = mc->PredictComponentSyst(calc,shift,
233  Flavors::kNuEToNuE,
234  Current::kCC,
235  Sign::kBoth).ToTHX(pot, force1D);
236  hAllBeamNue->SetNameTitle((hist_name+"_AllBeamNue").c_str(), (hist_name+"_AllBeamNue").c_str());
237  ret .push_back(hAllBeamNue);
238 
239  TH1* hBeamNue = mc->PredictComponentSyst(calc,shift,
240  Flavors::kNuEToNuE,
241  Current::kCC,
242  Sign::kNu).ToTHX(pot, force1D);
243  hBeamNue->SetNameTitle((hist_name+"_BeamNue").c_str(), (hist_name+"_BeamNue").c_str());
244  ret .push_back(hBeamNue);
245 
246  TH1* hBeamNuebar = mc->PredictComponentSyst(calc,shift,
247  Flavors::kNuEToNuE,
248  Current::kCC,
249  Sign::kAntiNu).ToTHX(pot, force1D);
250  hBeamNuebar->SetNameTitle((hist_name+"_BeamNuebar").c_str(), (hist_name+"_BeamNuebar").c_str());
251  ret .push_back(hBeamNuebar);
252 
253  TH1* hAllSignalNue = mc->PredictComponentSyst(calc,shift,
254  Flavors::kNuMuToNuE,
255  Current::kCC,
256  Sign::kBoth).ToTHX(pot, force1D);
257  hAllSignalNue->SetNameTitle((hist_name+"_AllSignalNue").c_str(), (hist_name+"_AllSignalNue").c_str());
258  ret .push_back(hAllSignalNue);
259 
260  TH1* hSignalNue = mc->PredictComponentSyst(calc,shift,
261  Flavors::kNuMuToNuE,
262  Current::kCC,
263  Sign::kNu).ToTHX(pot, force1D);
264  hSignalNue->SetNameTitle((hist_name+"_SignalNue").c_str(), (hist_name+"_SignalNue").c_str());
265  ret .push_back(hSignalNue);
266 
267  TH1* hSignalNuebar = mc->PredictComponentSyst(calc,shift,
268  Flavors::kNuMuToNuE,
269  Current::kCC,
270  Sign::kAntiNu).ToTHX(pot, force1D);
271  hSignalNuebar->SetNameTitle((hist_name+"_SignalNuebar").c_str(), (hist_name+"_SignalNuebar").c_str());
272  ret .push_back(hSignalNuebar);
273 
274  return ret;
275  }
276 
277  //----------------------------------------------------------------------
278  std::vector<TH1*> GetMCTotalForSystShifts(const IPrediction* mc,
279  osc::IOscCalc* calc,
280  const ISyst* syst,
281  std::string hist_base_name,
282  double pot,
283  bool force1D)
284  {
285  std::vector<TH1*> ret;
286  std::string syst_name = syst->ShortName();
287  for (int i = -3; i < 4; ++i){
288  SystShifts s(syst, double(i));
289  TH1* hTotal = mc->PredictSyst(calc, s).ToTHX(pot, force1D);
290  hTotal->SetNameTitle((hist_base_name+"_total_"+syst_name+"_"+std::to_string(i)).c_str(),
291  (hist_base_name+"_total_"+syst_name+"_"+std::to_string(i)).c_str());
292  ret .push_back(hTotal);
293  }
294  return ret;
295  }
296 
297 
298  //----------------------------------------------------------------------
299  void DataMCRatio(const Spectrum& data,
300  const IPrediction* mc,
301  osc::IOscCalc* calc,
302  double miny, double maxy)
303  {
304  DataMCRatio(data, mc->Predict(calc), miny, maxy);
305  }
306 
307  //----------------------------------------------------------------------
308  void DataMCRatio(const Spectrum& data,
309  const Spectrum& mc,
310  double miny, double maxy)
311  {
312  Ratio ratio(data, mc);
313  TH1* h = ratio.ToTH1();
314  h->GetYaxis()->SetTitle("Data / MC");
315  h->GetYaxis()->SetRangeUser(miny, maxy);
316  h->Draw();
317 
318  TLine* one = new TLine(h->GetXaxis()->GetXmin(), 1,
319  h->GetXaxis()->GetXmax(), 1);
320  one->SetLineWidth(2);
321  one->SetLineStyle(7);
322  one->Draw();
323 
324  gPad->Update();
325  }
326 
327  //----------------------------------------------------------------------
328  void DataMCAreaNormalizedRatio(const Spectrum& data,
329  const IPrediction* mc,
330  osc::IOscCalc* calc,
331  double miny, double maxy)
332  {
333  DataMCAreaNormalizedRatio(data, mc->Predict(calc), miny, maxy);
334  }
335 
336  //----------------------------------------------------------------------
337  void DataMCAreaNormalizedRatio(const Spectrum& data,
338  const Spectrum& mc,
339  double miny, double maxy)
340  {
341  Spectrum mcScaled = mc;
342  mcScaled.OverridePOT(mcScaled.POT() * mc.Integral(1) / data.Integral(1));
343 
344  DataMCRatio(data, mcScaled, miny, maxy);
345  }
346 
347  //----------------------------------------------------------------------
348  void RatioPlot(const Spectrum& data,
349  const Spectrum& expected,
350  const Spectrum& fit,
351  double miny, double maxy)
352  {
353  Ratio fitRatio(fit, expected);
354  Ratio dataRatio(data, expected);
355 
356  TH1* h = fitRatio.ToTH1();
357  h->GetYaxis()->SetTitle("Ratio to expectation");
358  h->GetYaxis()->SetRangeUser(miny, maxy);
359  h->SetLineColor(kTotalMCColor);
360  h->Draw("][");
361 
362  h = dataRatio.ToTH1();
363  h->SetMarkerStyle(kFullCircle);
364  h->Draw("ep same");
365 
366  TLine* one = new TLine(h->GetXaxis()->GetXmin(), 1,
367  h->GetXaxis()->GetXmax(), 1);
368  one->SetLineWidth(2);
369  one->SetLineStyle(7);
370  one->Draw();
371 
372  gPad->Update();
373  }
374 
375 
376  //----------------------------------------------------------------------
377  void PlotWithSystErrorBand(IPrediction* pred,
378  const std::vector<const ISyst*>& systs,
379  osc::IOscCalc* calc,
380  double pot,
381  int col, int errCol, float headroom,
382  bool newaxis)
383  {
384 
385  Spectrum nom = pred->Predict(calc);
386 
387  std::vector<Spectrum> ups, dns;
388 
389  for(const ISyst* syst: systs){
390  SystShifts shifts;
391  shifts.SetShift(syst, +1);
392  ups.push_back(pred->PredictSyst(calc, shifts));
393  shifts.SetShift(syst, -1);
394  dns.push_back(pred->PredictSyst(calc, shifts));
395  }
396 
397  PlotWithSystErrorBand(nom, ups, dns, pot, col, errCol, headroom, newaxis);
398 
399 
400  }
401 
402 
403  //----------------------------------------------------------------------
404  void PlotWithSystErrorBand(const Spectrum& nominal,
405  const std::vector<Spectrum>& upShifts,
406  const std::vector<Spectrum>& downShifts,
407  double pot, int col, int errCol,
408  float headroom, bool newaxis)
409  {
410  if(col == -1){
411  col = kTotalMCColor;
412  errCol = kTotalMCErrorBandColor;
413  }
414  else if(errCol == -1) errCol = col-7; // hopefully a lighter version
415 
416  TH1* nom = nominal.ToTH1(pot);
417 
418  std::vector<TH1*> ups, dns;
419 
420  for(const auto& upShift:upShifts) ups.push_back(upShift.ToTH1(pot));
421  for(const auto& downShift:downShifts) dns.push_back(downShift.ToTH1(pot));
422 
423  nom->SetLineColor(col);
424  nom->GetXaxis()->CenterTitle();
425  nom->GetYaxis()->CenterTitle();
426  if(newaxis) nom->Draw("hist ]["); // Set the axes up
427 
428  double yMax = nom->GetBinContent(nom->GetMaximumBin());
429 
430  TGraphAsymmErrors* g = new TGraphAsymmErrors;
431 
432  for(int binIdx = 0; binIdx < nom->GetNbinsX()+2; ++binIdx){
433  const double y = nom->GetBinContent(binIdx);
434  g->SetPoint(binIdx, nom->GetXaxis()->GetBinCenter(binIdx), y);
435 
436  const double w = nom->GetXaxis()->GetBinWidth(binIdx);
437 
438  double errUp = 0;
439  double errDn = 0;
440 
441  for(unsigned int systIdx = 0; systIdx < ups.size(); ++systIdx){
442  double hi = ups[systIdx]->GetBinContent(binIdx)-y;
443  double lo = y-dns[systIdx]->GetBinContent(binIdx);
444 
445  if(lo <= 0 && hi <= 0) std::swap(lo, hi);
446 
447  errUp += hi*hi;
448  errDn += lo*lo;
449 
450  // TODO: what happens if they're both high or both low?
451  } // end for systIdx
452 
453 
454  g->SetPointError(binIdx, w/2, w/2, sqrt(errDn), sqrt(errUp));
455  } // end for i
456 
457  g->SetFillColor(errCol);
458  g->Draw("e2 same");
459  g->GetYaxis()->SetRangeUser(0, headroom*yMax);
460  nom->GetYaxis()->SetRangeUser(0, headroom*yMax);
461 
462  nom->Draw("hist ][ same");
463 
464  for(TH1* up: ups) delete up;
465  for(TH1* dn: dns) delete dn;
466  }
467 
468  //----------------------------------------------------------------------
469  THStack* ToTHStack(const std::vector<std::pair<Spectrum, Color_t>>& s,
470  double pot)
471  {
472  THStack* ret = new THStack;
473  for(auto it: s){
474  TH1* h = it.first.ToTH1(pot);
475  h->SetFillColor(it.second);
476  ret->Add(h, "hist");
477  }
478  return ret;
479  }
480 
481  //----------------------------------------------------------------------
482  /// Helper for \ref AutoPlaceLegend
483  double PointDistanceToBox(double x, double y,
484  double x0, double y0, double x1, double y1)
485  {
486  // Inside
487  if(x > x0 && x < x1 && y > y0 && y < y1) return 0;
488 
489  // Corners
490  double d = util::sqr(x-x0)+util::sqr(y-y0);
491  d = std::min(d, util::sqr(x-x1)+util::sqr(y-y0));
492  d = std::min(d, util::sqr(x-x1)+util::sqr(y-y1));
493  d = std::min(d, util::sqr(x-x0)+util::sqr(y-y1));
494 
495  // Top and bottom edges
496  if(x > x0 && x < x1){
497  d = std::min(d, util::sqr(y-y0));
498  d = std::min(d, util::sqr(y-y1));
499  }
500  // Left and right
501  if(y > y0 && y < y1){
502  d = std::min(d, util::sqr(x-x0));
503  d = std::min(d, util::sqr(x-x1));
504  }
505 
506  return d;
507  }
508 
509  //----------------------------------------------------------------------
510  TLegend* AutoPlaceLegend(double dx, double dy, double yPin)
511  {
512  gPad->Update();
513 
514  // Convert requested width and height into physics coordinates
515  dx *= (gPad->GetX2()-gPad->GetX1());
516  dy *= (gPad->GetY2()-gPad->GetY1());
517 
518  // Range of axes in physics units
519  const double x0 = gPad->GetUxmin();
520  const double x1 = gPad->GetUxmax();
521  const double y0 = gPad->GetUymin();
522  const double y1 = gPad->GetUymax();
523 
524  const double X = x1-x0;
525  const double Y = y1-y0;
526 
527  double bestd = 0;
528  double bestx = 0;
529  double besty = 0;
530  // If we want to pin Y pos, set it to that now.
531  if(yPin >= 0) besty = yPin * (gPad->GetY2()-gPad->GetY1());
532 
533  for(bool fallback: {false, true}){
534  for(double x = x0+dx/2; x < x1-dx/2; x += X/50){
535  for(double y = y0+dy/2; y < y1-dy/2; y += Y/50){
536  double d = 999999;
537 
538  // Repel from edges
539  d = std::min(d, util::sqr((x-dx/2-x0)/X));
540  d = std::min(d, util::sqr((x+dx/2-x1)/X));
541  d = std::min(d, util::sqr((y-dy/2-y0)/Y));
542  d = std::min(d, util::sqr((y+dy/2-y1)/Y));
543  if(d < bestd) continue;
544 
545  TIter next(gPad->GetListOfPrimitives());
546  while(TObject* obj = next()){
547 
548  if(!obj->InheritsFrom(TH1::Class())) continue;
549  if( obj->InheritsFrom(TH2::Class())) continue;
550 
551  TH1* h = (TH1*)obj;
552 
553  for(int n = 1; n <= h->GetNbinsX(); ++n){
554  const double px = h->GetBinCenter(n);
555  const double py = h->GetBinContent(n);
556 
557  if(fallback){
558  d = std::min(d, util::sqr((px-x)/X)+util::sqr((py-y)/Y));
559  }
560  else{
561  d = std::min(d, PointDistanceToBox(px/X, py/Y,
562  (x-dx/2)/X, (y-dy/2)/Y,
563  (x+dx/2)/X, (y+dy/2)/Y));
564  }
565  if(d < bestd) break;
566  }
567  if(d < bestd) break;
568  } // end while
569 
570  if(d > bestd){
571  bestd = d;
572  bestx = x;
573  // Update Y if we're not pinning it.
574  if (yPin < 0) besty = y;
575  }
576  } // end for y
577  } // end for x
578 
579  if(bestd != 0) break; // If we always collide, have to do fallback
580  } // end for fallback
581 
582  // Convert to pad coordinates
583  const double nx = (bestx-gPad->GetX1())/(gPad->GetX2()-gPad->GetX1());
584  const double ny = (besty-gPad->GetY1())/(gPad->GetY2()-gPad->GetY1());
585 
586  const double ndx = dx/(gPad->GetX2()-gPad->GetX1());
587  const double ndy = dy/(gPad->GetY2()-gPad->GetY1());
588 
589  return new TLegend(nx-ndx/2, ny-ndy/2, nx+ndx/2, ny+ndy/2);
590  }
591 
592  //----------------------------------------------------------------------
593  /*
594  void CountingExperimentErrorBarChart(const std::map<std::string, double>& systbars, double statErr, bool bkgdOrSig, bool shortchart)
595  {
596  std::map<std::string, double> systs = (shortchart) ? SumNueSecondAnaSysts(systbars) : systbars;
597 
598  int nbins = systs.size()+1;
599  // Systematics are on top, total systematic and (optional) statistical
600  // error bar are at the bottom.
601  int firstSystBin = 1;
602  if(statErr){
603  ++nbins;
604  ++firstSystBin;
605  }
606 
607  // Sum all the systematics in qaudrature
608  double systTot = 0;
609  for(auto it: systs) systTot += it.second*it.second;
610  systTot = sqrt(systTot);
611 
612  // x-axis range is a multiple of 5% with at least 2% padding
613  double xrange = 0;
614  while(xrange < std::max(systTot, statErr)+2) xrange += 5;
615 
616  TH2* axes = new TH2F("", ";Background uncertainty (%)",
617  10, -xrange, +xrange, nbins, 0, nbins);
618  axes->GetXaxis()->CenterTitle();
619  if(bkgdOrSig) axes->GetXaxis()->SetTitle("Signal uncertainty (%)");
620 
621  // Put the systematics in a vector so we can sort them by size
622  std::vector<std::pair<double, std::string>> systList;
623  for(auto it: systs) systList.push_back(std::make_pair(it.second, it.first));
624  std::sort(systList.begin(), systList.end());
625 
626  // Label the y axis
627  TAxis* yax = axes->GetYaxis();
628  for(unsigned int i = 0; i < systList.size(); ++i)
629  yax->SetBinLabel(i+firstSystBin+1, systList[i].second.c_str());
630 
631  yax->SetBinLabel(firstSystBin, "Total syst. error");
632  if(statErr) yax->SetBinLabel(1, "Statistical error");
633 
634  // Make the labels approximately match the x axis title
635  yax->SetLabelSize(.06);
636 
637  axes->Draw();
638 
639  // Blue boxes for each systematic
640  for(unsigned int i = 0; i < systList.size(); ++i){
641  TBox* box = new TBox(-systList[i].first, i+firstSystBin+.2,
642  +systList[i].first, i+firstSystBin+.8);
643  std::cout<<"ERROR = "<<systList[i].first<<std::endl;
644  box->SetFillColor(kAzure-8);
645  box->Draw();
646  }
647 
648  // Total systematic error
649  TBox* syst = new TBox(-systTot, firstSystBin-.8,
650  +systTot, firstSystBin-.2);
651  std::cout<<"TOTAL SYSTEMATIC ERROR = "<<systTot<<std::endl;
652  syst->SetFillColor(kAzure-8);
653  syst->Draw();
654 
655  // Red box for statistical error
656  if(statErr){
657  TBox* stat = new TBox(-statErr, .2, +statErr, .8);
658  stat->SetFillColor(kRed-7);
659  stat->Draw();
660  }
661 
662  // Separate individual systematics from the total and statistical errors
663  TLine* div = new TLine(-xrange, firstSystBin, +xrange, firstSystBin);
664  div->SetLineWidth(2);
665  div->Draw();
666 
667  // Vertical line marking the symmetry point
668  TLine* zero = new TLine(0, 0, 0, nbins);
669  zero->SetLineStyle(7);
670  zero->SetLineWidth(2);
671  zero->Draw();
672 
673  // Leave enough space for the systematic labels
674  gPad->SetLeftMargin(.25);
675  }
676  */
677  //----------------------------------------------------------------------
678  TGraphAsymmErrors* GraphWithPoissonErrors(const TH1* h, bool noErrorsXaxis, bool drawEmptyBins)
679  {
680  TGraphAsymmErrors* gr = new TGraphAsymmErrors(h);
681 
682  TFeldmanCousins fc(0.6827);//1 sigma
683 
684  for(int i = 0; i < h->GetNbinsX(); ++i){
685  double x, y;
686  gr->GetPoint(i, x, y);
687 
688  if ( drawEmptyBins || y!=0 ){
689  if ( y < 50 ) gr->SetPointEYlow(i, y-fc.CalculateLowerLimit(y,0));
690  else gr->SetPointEYlow(i,sqrt(y));
691  if ( y < 30 ) gr->SetPointEYhigh(i,fc.CalculateUpperLimit(y,0)-y);
692  else gr->SetPointEYhigh(i,sqrt(y));
693  }
694 
695  if(noErrorsXaxis){
696  gr->SetPointEXlow(i, 0);
697  gr->SetPointEXhigh(i, 0);
698  } // Do not use bin width as X-error for points
699  }
700 
701  return gr;
702  }
703  //----------------------------------------------------------------------
704  TGraph* ShadeBetweenHistograms(TH1* hmin, TH1* hmax)
705  {
706  int n = hmin->GetNbinsX();
707  TGraph* gr = new TGraph(4*n);
708 
709  for(int i=1; i<=n; i++)
710  {
711  double xdown = hmax->GetBinLowEdge(i);
712  double xup = hmax->GetBinLowEdge(i+1);
713  double ydown = hmin->GetBinContent(i);
714  double yup = hmax->GetBinContent(i);
715 
716  gr->SetPoint(2*(i-1), xdown, ydown);
717  gr->SetPoint(2*(i-1)+1, xup, ydown);
718  gr->SetPoint(4*n-2*i, xup, yup);
719  gr->SetPoint(4*n-2*i+1, xdown, yup);
720  }
721 
722  gr->SetFillColor(hmin->GetLineColor() - 9); // In principle this is lighter
723  gr->SetLineColor(hmin->GetLineColor() - 9); // In case one does Draw("lf")
724 
725  return gr;
726  }
727  //----------------------------------------------------------------------
728  TGraphAsymmErrors * ProfileQuantile(const TH2 * hist,
729  const std::string & axisName,
730  const std::string & graphName,
731  const std::pair<double, double> & quantileDivisions)
732  {
733  if (hist->GetDimension() != 2)
734  throw std::runtime_error(Form("Can't profile a histogram with other than 2 dimensions. Yours is a %d-D histogram...", hist->GetDimension()));
735 
736  const TAxis * axis = nullptr;
737  std::function<TH1D*(const char *, Int_t, Int_t, Option_t*)> projectionMethod;
738  if (axisName == "x" || axisName == "X")
739  {
740  axis = hist->GetXaxis();
741  projectionMethod = std::bind(&TH2::ProjectionY, hist, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4);
742  }
743  else if (axisName == "y" || axisName == "Y")
744  {
745  axis = hist->GetYaxis();
746  projectionMethod = std::bind(&TH2::ProjectionX, hist, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4);
747  }
748  else
749  throw std::runtime_error( Form("Can't profile onto unknown axis: '%s'", axisName.c_str()) );
750 
751 
752  std::vector<double> points_x, points_y, errors_x_low, errors_x_high, errors_y_low, errors_y_high;
753 
754  double quantiles[2] = {0, 0};
755  double _quantileDivisions[2] = {quantileDivisions.first, quantileDivisions.second};
756  for (int binNum = 1; binNum <= axis->GetNbins(); binNum++) // remember, 0 is underflow and Nbins+1 is overflow...
757  {
758  std::unique_ptr<TH1D> proj ( projectionMethod(Form("%s_tmp_%d", hist->GetName(), int(std::time(nullptr))), binNum, binNum, "") ); // randomish name, just to be safe
759 
760  double y = proj->GetMean();
761  points_x.push_back(axis->GetBinCenter(binNum));
762  points_y.push_back(y);
763 
764  // for now, make the errors for an empty projection 0
765  unsigned int n_qs = 0;
766  if (proj->Integral() == 0)
767  {
768  n_qs = 2;
769  quantiles[0] = quantiles[1] = 0;
770  }
771  else
772  n_qs = proj->GetQuantiles(2, quantiles, _quantileDivisions);
773  if (n_qs != 2)
774  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) );
775 
776  double binWidth = axis->GetBinWidth(binNum);
777  errors_x_low.push_back(binWidth/2.);
778  errors_x_high.push_back(binWidth/2.);
779  errors_y_low.push_back(y-quantiles[0]);
780  errors_y_high.push_back(quantiles[1]-y);
781  }
782 
783  TGraphAsymmErrors * outGraph = new TGraphAsymmErrors(
784  points_x.size(),
785  &points_x[0],
786  &points_y[0],
787  &errors_x_low[0],
788  &errors_x_high[0],
789  &errors_y_low[0],
790  &errors_y_high[0]
791  );
792  std::string name = (graphName.size()) ? graphName : Form("%s_quantile_errors", hist->GetName());
793  outGraph->SetName(name.c_str());
794 
795  return outGraph;
796 
797  }
798 } // namespace
ana::ToTHStack
THStack * ToTHStack(const std::vector< std::pair< Spectrum, Color_t >> &s, double pot)
Can call like ToTHStack({{h1, kRed}, {h2, kBlue}}, pot)
Definition: Plots.cxx:469
ana::DataMCRatio
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:299
ana::Spectrum::ToTH1
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.
ana::SystShifts::SetShift
void SetShift(const ISyst *syst, double shift)
Shifts are 0=nominal, -1,+1 = 1 sigma shifts.
Definition: SystShifts.cxx:47
ana::Sign::kAntiNu
Antineutrinos-only.
Definition: IPrediction.h:49
ana::Ratio::ToTH1
TH1D * ToTH1(Color_t col=kBlack, Style_t style=kSolid) const
Definition: Ratio.cxx:116
x
process_name opflash particleana ie x
Definition: run_opflash_electron.fcl:90
ana::Flavors::kNuMuToNuE
(&#39; appearance&#39;)
Definition: IPrediction.h:17
ana::DataMCComparison
TH1 * DataMCComparison(const Spectrum &data, const Spectrum &mc)
Definition: Plots.cxx:32
ana::PlotWithSystErrorBand
void PlotWithSystErrorBand(IPrediction *pred, const std::vector< const ISyst * > &systs, osc::IOscCalc *calc, double pot, int col, int errCol, float headroom, bool newaxis)
Plot prediction with +/-1sigma error band.
Definition: Plots.cxx:377
ana::SystShifts
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:16
ana::Spectrum::OverridePOT
void OverridePOT(double newpot)
DO NOT USE UNLESS YOU ARE 110% CERTAIN THERE ISN&#39;T A BETTER WAY!
Definition: Spectrum.h:295
ana::Flavors::kNuEToNuE
(&#39;beam &#39;)
Definition: IPrediction.h:14
osc::_IOscCalc
Definition: Plots.h:17
ana::IPrediction::PredictComponentSyst
virtual Spectrum PredictComponentSyst(osc::IOscCalc *calc, const SystShifts &syst, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const
Definition: IPrediction.cxx:72
ana::GraphWithPoissonErrors
TGraphAsymmErrors * GraphWithPoissonErrors(const TH1 *h, bool noErrorsXaxis, bool drawEmptyBins)
Calculate statistical errors appropriate for small Poisson numbers.
Definition: Plots.cxx:678
ana
process_name opflashCryoW ana
Definition: stage0_icarus_light_only.fcl:60
ana::Spectrum::Integral
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
X
then echo echo For and will not be changed by echo further linking echo echo B echo The symbol is in the uninitialized data multiple common symbols may appear with the echo same name If the symbol is defined the common echo symbols are treated as undefined references For more echo details on common see the discussion of warn common echo in *Note Linker see the discussion of warn common echo in *Note Linker such as a global int variable echo as opposed to a large global array echo echo I echo The symbol is an indirect reference to another symbol This echo is a GNU extension to the a out object file format which is echo rarely used echo echo N echo The symbol is a debugging symbol echo echo R echo The symbol is in a read only data section echo echo S echo The symbol is in an uninitialized data section for small echo objects echo echo T echo The symbol is in the the normal defined echo symbol is used with no error When a weak undefined symbol echo is linked and the symbol is not the value of the echo weak symbol becomes zero with no error echo echo W echo The symbol is a weak symbol that has not been specifically echo tagged as a weak object symbol When a weak defined symbol echo is linked with a normal defined the normal defined echo symbol is used with no error When a weak undefined symbol echo is linked and the symbol is not the value of the echo weak symbol becomes zero with no error echo echo echo The symbol is a stabs symbol in an a out object file In echo this the next values printed are the stabs other echo the stabs desc and the stab type Stabs symbols are echo used to hold debugging information For more echo see *Note or object file format specific echo echo For Mac OS X
Definition: find_global_symbol.sh:107
ana::kTotalMCErrorBandColor
const Color_t kTotalMCErrorBandColor
Definition: sbnana/sbnana/CAFAna/Analysis/Style.h:14
g
BEGIN_PROLOG g
Definition: decoderTools_icarus.fcl:7
shift
shift
Definition: fcl_checks.sh:26
util::sqr
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
ana::ISyst
Encapsulate code to systematically shift a caf::StandardRecord.
Definition: ISyst.h:14
ana::kNumuBackgroundColor
const Color_t kNumuBackgroundColor
Definition: sbnana/sbnana/CAFAna/Analysis/Style.h:22
ana::ISyst::ShortName
virtual std::string ShortName() const final
The name printed out to the screen.
Definition: ISyst.h:30
ana::ProjectionY
void ProjectionY(TH2D *from, TH1D *to)
Helper for WeightingVariable.
Definition: ReweightableSpectrum.cxx:294
ana::ProfileQuantile
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:728
ana::Spectrum
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:30
ana::GetMCTotalForSystShifts
std::vector< TH1 * > GetMCTotalForSystShifts(const IPrediction *mc, osc::IOscCalc *calc, const ISyst *syst, std::string hist_base_name, double pot, bool force1D)
Definition: Plots.cxx:278
ana::IPrediction::Predict
virtual Spectrum Predict(osc::IOscCalc *calc) const =0
h
while getopts h
Definition: getStatistics.sh:13
lar::dump::vector
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:265
ana::Current::kCC
Charged-current interactions.
Definition: IPrediction.h:38
y
process_name opflash particleana ie ie y
Definition: run_opflash_electron.fcl:91
ana::RatioPlot
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:348
ana::DataMCComparisonAreaNormalized
TH1 * DataMCComparisonAreaNormalized(const Spectrum &data, const Spectrum &mc)
Definition: Plots.cxx:62
ana::kBeamNueBackgroundColor
const Color_t kBeamNueBackgroundColor
Definition: sbnana/sbnana/CAFAna/Analysis/Style.h:18
ana::ShadeBetweenHistograms
TGraph * ShadeBetweenHistograms(TH1 *hmin, TH1 *hmax)
Definition: Plots.cxx:704
ana::Spectrum::ToTHX
TH1 * ToTHX(double exposure, bool force1D=false, EExposureType expotype=kPOT) const
Function decides what is the appropriate projection based on fBins, and does that.
ana::GetMCTotal
TH1 * GetMCTotal(const IPrediction *mc, osc::IOscCalc *calc, std::string hist_name, double pot, bool force1D)
Definition: Plots.cxx:158
ana::IPrediction::PredictSyst
virtual Spectrum PredictSyst(osc::IOscCalc *calc, const SystShifts &syst) const
Definition: IPrediction.cxx:62
ana::ProjectionX
void ProjectionX(TH2D *from, TH1D *to)
Helper for Unweighted.
Definition: ReweightableSpectrum.cxx:255
ana::GetMCSystComponents
std::vector< TH1 * > GetMCSystComponents(const IPrediction *mc, osc::IOscCalc *calc, const SystShifts &shift, std::string hist_name, double pot, bool force1D)
Definition: Plots.cxx:176
ana::Spectrum::POT
double POT() const
Definition: Spectrum.h:289
ana::Ratio
Represent the ratio between two spectra.
Definition: Ratio.h:8
gen_crt_frags.w
tuple w
Definition: gen_crt_frags.py:1416
ana::Sign::kNu
Neutrinos-only.
Definition: IPrediction.h:48
ana::kNoShift
const SystShifts kNoShift
Definition: SystShifts.h:61
ana::PointDistanceToBox
double PointDistanceToBox(double x, double y, double x0, double y0, double x1, double y1)
Helper for AutoPlaceLegend.
Definition: Plots.cxx:483
one
then echo Cowardly refusing to create a new FHiCL file with the same name as the original one('${SourceName}')." >&2 exit 1 fi echo "'$
Definition: makeNoOverburdenVersionOf.sh:65
ana::DataMCAreaNormalizedRatio
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:328
px
BEGIN_PROLOG px
Definition: filemuons.fcl:10
ana::DataMCComparisonComponents
TH1 * DataMCComparisonComponents(const Spectrum &data, const IPrediction *mc, osc::IOscCalc *calc)
Plot MC broken down into flavour components, overlayed with data.
Definition: Plots.cxx:94
py
process_name physics producers generator physics producers generator physics producers generator py
Definition: prod_muon_2GeV_d45_Collection.fcl:21
icarus::trigger::to_string
std::string to_string(WindowPattern const &pattern)
Definition: WindowPattern.h:115
s
then echo File list $list not found else cat $list while read file do echo $file sed s
Definition: file_to_url.sh:60
ana::GetMCComponents
std::vector< TH1 * > GetMCComponents(const IPrediction *mc, osc::IOscCalc *calc, std::string hist_name, double pot, bool force1D)
Definition: Plots.cxx:167
ana::GetMCSystTotal
TH1 * GetMCSystTotal(const IPrediction *mc, osc::IOscCalc *calc, const SystShifts &shift, std::string hist_name, double pot, bool force1D)
Definition: Plots.cxx:146
ana::Current::kNC
Neutral-current interactions.
Definition: IPrediction.h:39
ana::IPrediction::PredictComponent
virtual Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
ana::Sign::kBoth
Both neutrinos and antineutrinos.
Definition: IPrediction.h:51
ana::IPrediction
Standard interface to all prediction techniques.
Definition: IPrediction.h:58
name
then echo fcl name
Definition: sbndpoms_genfclwithrunnumber_maker.sh:220
ana::AutoPlaceLegend
TLegend * AutoPlaceLegend(double dx, double dy, double yPin)
Create a legend, maximizing distance from all histograms.
Definition: Plots.cxx:510
ana::kTotalMCColor
const Color_t kTotalMCColor
Definition: sbnana/sbnana/CAFAna/Analysis/Style.h:13
ana::Flavors::kAll
All neutrinos, any flavor.
Definition: IPrediction.h:25
ana::kNCBackgroundColor
const Color_t kNCBackgroundColor
Definition: sbnana/sbnana/CAFAna/Analysis/Style.h:17
Generated by   doxygen 1.8.5