systematics_summary_utils.h
Go to the documentation of this file.
4 
5 #include "OscLib/IOscCalc.h"
6 #include "OscLib/OscCalcDumb.h"
7 
8 //#include "3FlavorAna/Ana2020/joint_fit_style_tools.h"
11 
13 
14 #include "TCanvas.h"
15 #include "TGraphAsymmErrors.h"
16 #include "TPaveText.h"
17 #include "TFeldmanCousins.h"
18 
19 #include <algorithm>
20 #include <string>
21 #include <iostream>
22 #include <iomanip>
23 
24 TString FixSystName(TString mystring){
25  //There must be a better way
26  std::vector<TString> in = {"#",".","_"};
27  std::vector<TString> out = {" "," "," "};
28  for(unsigned int i=0;i<in.size();i++)
29  mystring.ReplaceAll(in[i],out[i]);
30  return mystring;
31 }
32 
33 // Summary tables' columns definitios
34 struct ColumnDef
35 {
40 };
41 
42 // Systematic categories for Ana2020
52 };
53 
54 // Names of the systematic categories
56  {"Near-Far Uncor.",
57  "Neutrino Cross Sections",
58  "Detector Calibration",
59  "Beam Flux",
60  "Detector Response",
61  "Neutron Uncertainty",
62  "Lepton Reconstruction"
63 };
64 
65 // Line colors for syst categories for ratioband plots
67  kGreen+1, // Near-Far
68  kAzure-3, // Xsecs
69  kMagenta-3, // Calib
70  0, // Beam, this syst is so low there's no point in plotting it
71  kOrange+8, // Det Response
72  kGray+2, // Neutron
73  kCyan-7 // Lepton Reco
74 };
75 
77  // ---- Print out my info to start the document...
78  std::cout << "\n The latex code follows. Note that you need to do some post editing;"
79  << "\n\t Remove any cases of #1, #2, etc"
80  << "\n\t Change cases of #bar{#nu_{e}}, and all sub cases"
81  << "\n\t Change cases of |#vec{q}|"
82  << "\n\t Change cases of Q^2"
83  << std::endl;
84 
85  std::cout<<"\n\n" << std::endl;
86  std::cout<<"\\documentclass[10pt]{article}"<<std::endl;
87  std::cout<<"\\usepackage[usenames]{color} %used for font color"<<std::endl;
88  std::cout<<"\\usepackage{amssymb} %maths"<<std::endl;
89  std::cout<<"\\usepackage{amsmath} %maths"<<std::endl;
90  std::cout<<"\\usepackage[utf8]{inputenc} %useful to type directly diacritic characters"<<std::endl;
91  std::cout<<"\\usepackage{color, colortbl}"<<std::endl;
92  std::cout<<"\\usepackage{graphicx}"<<std::endl;
93  std::cout<<"\\usepackage[dvipsnames]{xcolor}"<<std::endl;
94  std::cout<<"\\usepackage{pdflscape}"<<std::endl;
95 
96  std::cout<<"\\usepackage{booktabs}"<<std::endl;
97 
98  std::cout<<"\\title{Systematics Summary for $\\nu_\\mu$ FHC 2020}" << std::endl;
99  std::cout<<"\\author{T. Nosek}"<<std::endl;
100  std::cout<<"\\date{\\today}"<<std::endl;
101  std::cout<<"\\setlength{\\textheight}{20.5cm}"<<std::endl;
102 
103  std::cout<<"\\begin{document}"<<std::endl;
104  std::cout<<"\\maketitle"<<std::endl;
105  std::cout<<"\\vfill"<<std::endl;
106  std::cout<<"\\tableofcontents"<<std::endl;
107  std::cout<<"\\centering"<<std::endl;
108 
109  //std::cout<<"\n\\begin{document}"<<std::endl;
110  //if (!IsNuMu) std::cout<<"\\begin{landscape}"<<std::endl;
111 }
112 
113 void drawLabel(std::string s,double x=0.3,double y=0.93){
114  TLatex *tex = new TLatex(x,y,s.c_str());
115  tex->SetTextSize(1./30);
116  tex->SetNDC();
117  tex->SetTextAlign(12);
118  tex->Draw();
119 }
120 
121 void drawBeamLabel(bool isFHC, double x=0.35){
122  std::string s;
123  if (isFHC) s = "#nu-beam";
124  else s = "#bar{#nu}-beam";
125  TLatex *tex = new TLatex(x, 0.93, s.c_str());
126  tex->SetTextColor(kGray+2);
127  tex->SetNDC();
128  tex->SetTextSize(2/40.);
129  tex->SetTextAlign(12);
130  tex->Draw();
131 }
132 
133 
134 void drawSampleLabel(std::string s, double x=.915, double y=.9){
135  TLatex *tex = new TLatex(x,y,s.c_str());
136  tex->SetTextSize(2./40.);
137  tex->SetNDC();
138  tex->SetTextAngle(270);
139  //tex->SetTextAlign(12);
140  tex->Draw();
141 }
142 
143 // Spectrum integral
144 double Integral(const Spectrum& s, double exp, double* err = 0, EExposureType expotype = kPOT)
145 {
146  return s.Integral(exp, err, expotype);
147 }
148 
150 {
151  TLatex* prelim = new TLatex(.9, .93, "NOvA Preliminary");
152  prelim->SetTextColor(kAzure);
153  prelim->SetNDC();
154  prelim->SetTextSize(2/40.);
155  prelim->SetTextAlign(32);
156  prelim->Draw();
157 }
158 
159 // Error band syst plots
160 void MyPlotWithSystErrorBand(TH1*& nom,
161  std::vector<TH1*>& ups,
162  std::vector<TH1*>& dns,
163  int col, int errCol,
164  float linewidth, bool newaxis)
165 {
166  if(col == -1){
167  col = kTotalMCColor;
168  errCol = kTotalMCErrorBandColor;
169  }
170 
171  else if(errCol == -1) errCol = col-7; // hopefully a lighter version
172  nom->SetLineColor(col);
173  nom->GetXaxis()->CenterTitle();
174  nom->GetYaxis()->CenterTitle();
175  nom->GetYaxis()->SetTitle("1#sigma Syst. Band Ratio");
176  if(newaxis) nom->Draw("hist ]["); // Set the axes up
177  double yMax = nom->GetBinContent(nom->GetMaximumBin());
178  TH1* hUp = (TH1*)nom->Clone();
179  TH1* hDn = (TH1*)nom->Clone();
180  for(int binIdx = 0; binIdx < nom->GetNbinsX()+2; ++binIdx){
181  const double y = nom->GetBinContent(binIdx);
182  if(y==0) nom->SetBinContent(binIdx,1);
183  double errUp = 0;
184  double errDn = 0;
185  for(unsigned int systIdx = 0; systIdx < ups.size(); ++systIdx){
186  double hi = ups[systIdx]->GetBinContent(binIdx)-y;
187  double lo = dns[systIdx]->GetBinContent(binIdx)-y;
188  // If both shifts are in the same direction use the larger
189  double min = std::min(hi,lo);
190  double max = std::max(hi,lo);
191  if(max < 0) max=0;
192  if(min > 0) min=0;
193  errUp += max*max;
194  errDn += min*min;
195  } // end for systIdx
196  hUp->SetBinContent(binIdx, y+sqrt(errUp));
197  if((y+sqrt(errUp)) > 1.5 || (y+sqrt(errUp)) < 1) hUp->SetBinContent(binIdx, 1);
198  hDn->SetBinContent(binIdx, y-sqrt(errDn));
199  if((y-sqrt(errDn)) < .1 || (y-sqrt(errDn)) > 1) hDn->SetBinContent(binIdx, 1);
200  } // end for i
201  hUp->SetLineColor(errCol);
202  hUp->SetLineWidth(linewidth);
203  hDn->SetLineColor(errCol);
204  hDn->SetLineWidth(linewidth);
205  hUp->Draw("hist ][ same");
206  hDn->Draw("hist ][ same");
207  // Use the yMax value or set the range by hand
208  //g->GetYaxis()->SetRangeUser(0.00001, headroom*yMax);
209  //g->GetYaxis()->SetRangeUser(0.8, 1.2);
210  //nom->GetYaxis()->SetRangeUser(0.0001, headroom*yMax);
211  nom->GetYaxis()->SetRangeUser(0.65, 1.35);
212  nom->Draw("hist ][ same");
213 }
214 
215 // bar chart to compare more predictions (like etrap vs noextrap)
217  std::vector<std::string> ylabels,
218  std::vector<std::vector<std::pair<double,double>>> errs,
219  std::vector<std::pair<double,double>> sum,
220  std::vector<std::string> legend,
221  std::vector<int> colors)
222 {
223  // add more colors as needed
224  if(colors.back()==-1) colors = {kAzure-8,kRed-7,kGreen+2};
225  else if(colors.back()==-2) colors = {kRed-4,kViolet-5,kViolet+10};
226 
227  int nbins = ylabels.size() + 1;
228  int firstBin = 1;
229 
230  double xrange = 0;
231  double max = 0;
232  for(unsigned int i = 0; i < sum.size(); ++i)
233  if(std::max(sum[i].first,sum[i].second) > max)
234  max = std::max(sum[i].first,sum[i].second);
235  while(xrange < max + 2 ) xrange += 5;
236 
237  TH2* axes = new TH2F("", (";" + xlabel).c_str(),
238  10, -xrange, +xrange, nbins, 0, nbins);
239 
240  axes->GetXaxis()->CenterTitle();
241 
242  TAxis* yax = axes->GetYaxis();
243  for(unsigned int i = 0; i < ylabels.size(); ++i)
244  yax->SetBinLabel(i+2, ylabels[i].c_str());
245  yax->SetBinLabel(1, "Systematic Uncertainty");
246  yax->SetLabelSize(.07);
247  //yax->SetLabelSize(.02);
248 
249  axes->Draw();
250 
251  TLegend *l = new TLegend(0.354,.894-.048*sum.size(),0.601,0.894);
252  l->SetFillStyle(0);
253 
254  for(unsigned int i = 0; i < errs.size(); ++i){
255  double size = .8/errs.size();
256  for(unsigned int j = 0; j < errs[i].size(); ++j){
257  TBox* box = new TBox(-errs[i][j].first, j+firstBin+.1+(errs.size()-i-1)*size,
258  +errs[i][j].second, j+firstBin+.1+(errs.size()-i)*size);
259  box->SetFillColor(colors[i]);
260  box->SetLineColor(colors[i]);
261  box->Draw();
262  }
263 
264  TBox* syst = new TBox(-sum[i].first, 0.1+(errs.size()-i-1)*size,
265  +sum[i].second, 0.1+(errs.size()-i)*size);
266  syst->SetFillColor(colors[i]);
267  syst->SetLineColor(colors[i]);
268  syst->Draw();
269 
270  l->AddEntry(syst,legend[i].c_str(),"f");
271  }
272 
273  l->Draw();
274 
275  TLine* div = new TLine(-xrange, 1, +xrange, 1);
276  div->SetLineWidth(2);
277  div->Draw();
278 
279  TLine* zero = new TLine(0, 0, 0, nbins);
280  zero->SetLineStyle(7);
281  zero->SetLineWidth(2);
282  zero->Draw();
283 
284  gPad->SetLeftMargin(.3);
285 }
286 
287 std::map<std::string, std::pair<double,double>> SumSysts(std::map<std::string, std::pair<double,double>> bars)
288 {
289  std::map<std::string, std::pair<double,double>> summedbars;
290  std::pair<double,double> summedsyst[kNumSystTypes];
291  for(int i = 0; i < kNumSystTypes; ++i) summedsyst[i] = {0,0};
292 
293  for(auto it: bars){
294  // Near-Far Uncor. systs
295  if( it.first.find("accept") != std::string::npos ||
296  it.first.find("michel") != std::string::npos ||
297  it.first.find("NormHorn") != std::string::npos ||
298  it.first.find("NormFHC") != std::string::npos ||
299  it.first.find("NormRHC") != std::string::npos ||
300  it.first.find("RockScale") != std::string::npos ||
301  it.first.find("cosmicScale") != std::string::npos )
302  {
303  //std::cerr << "NF diffs adding: -->" << it.first << std::endl;
304  summedsyst[kNFUncorSys].first += it.second.first*it.second.first;
305  summedsyst[kNFUncorSys].second+= it.second.second*it.second.second;
306  continue;
307  }
308 
309  // Lepton Reco systs
310  else if( it.first.find("MuE") != std::string::npos ||
311  it.first.find("LeptonAngle") != std::string::npos )
312  {
313  //std::cerr << "MuE scale adding: -->" << it.first << std::endl;
314  summedsyst[kLeptonRecoSys].first += it.second.first*it.second.first;
315  summedsyst[kLeptonRecoSys].second+= it.second.second*it.second.second;
316  continue;
317  }
318 
319  // Calibration systs
320  else if( it.first.find("Calib") != std::string::npos )
321  {
322  //std::cerr << "Detector calibration adding: -->" << it.first << std::endl;
323  summedsyst[kCalibSys].first += it.second.first*it.second.first;
324  summedsyst[kCalibSys].second+= it.second.second*it.second.second;
325  continue;
326  }
327 
328  // Neutron visible E systs
329  else if( it.first.find("NeutronEvis") != std::string::npos )
330  {
331  //std::cerr << "Neutron uncert adding: -->" << it.first << std::endl;
332  summedsyst[kNeutronSys].first += it.second.first*it.second.first;
333  summedsyst[kNeutronSys].second += it.second.second*it.second.second;
334  continue;
335  }
336 
337  // Detector response systs
338  else if( it.first.find("Cherenkov") != std::string::npos ||
339  it.first.find("Light Level") != std::string::npos )
340  {
341  //std::cerr << "Detector response adding: -->" << it.first << std::endl;
342  summedsyst[kBirksSys].first += it.second.first*it.second.first;
343  summedsyst[kBirksSys].second+= it.second.second*it.second.second;
344  continue;
345  }
346 
347  // Beam flux systs
348  else if( it.first.find("ppfx") != std::string::npos )
349  {
350  //std::cerr << "Beam flux adding: -->" << it.first << std::endl;
351  summedsyst[kBeamSys].first += it.second.first*it.second.first;
352  summedsyst[kBeamSys].second+= it.second.second*it.second.second;
353  continue;
354  }
355 
356  // Else Genie Xsec systs
357  else
358  {
359  //std::cerr << "Neutrino cross sections adding: -->" << it.first << std::endl;
360  summedsyst[kGenieSys].first += it.second.first*it.second.first;
361  summedsyst[kGenieSys].second+= it.second.second*it.second.second;
362  continue;
363  }
364 
365  std::cout<<"it.first == \""<<it.first<<"\" || "<<std::endl;
366  }
367 
368  for (int i = 0; i < kNumSystTypes; ++i){
369  summedbars[systTypes[i]] = {sqrt(summedsyst[i].first),sqrt(summedsyst[i].second)};
370  }
371 
372  return summedbars;
373 }
374 
375 std::vector<TH1*> MakeRatios(TH1* nom, std::vector<TH1*> shifts)
376 {
377  std::vector<TH1*> ret;
378  for(auto shift: shifts)
379  {
380  TH1* hTemp = (TH1*)shift->Clone();
381  hTemp->Divide(nom);
382  ret.push_back(hTemp);
383  }
384 
385  return ret;
386 }
387 
388 std::map<std::string, std::vector<TH1*>> SortSystHists(std::map<std::string, TH1*> shifts)
389 {
390  std::map<std::string, std::vector<TH1*> > ret;
391  std::vector<TH1*> tempRet[kNumSystTypes];
392 
393  for(auto shift: shifts)
394  {
395  if( shift.first.find("accept") != std::string::npos ||
396  shift.first.find("michel") != std::string::npos ||
397  shift.first.find("NormHorn") != std::string::npos ||
398  shift.first.find("NormFHC") != std::string::npos ||
399  shift.first.find("NormRHC") != std::string::npos ||
400  shift.first.find("RockScale") != std::string::npos ||
401  shift.first.find("cosmicScale") != std::string::npos )
402  {
403  std::cerr << "NF Uncor adding: -->" << shift.first << std::endl;
404  tempRet[kNFUncorSys].push_back(shift.second);
405  continue;
406  }
407 
408  // Lepton Reco systs
409  else if( shift.first.find("MuE") != std::string::npos ||
410  shift.first.find("LeptonAngle") != std::string::npos )
411  {
412  std::cerr << "MuE scale adding: -->" << shift.first << std::endl;
413  tempRet[kLeptonRecoSys].push_back(shift.second);
414  continue;
415  }
416 
417  // Calibration systs
418  else if( shift.first.find("Calib") != std::string::npos )
419  {
420  std::cerr << "Detector calibration adding: -->" << shift.first << std::endl;
421  tempRet[kCalibSys].push_back(shift.second);
422  continue;
423  }
424 
425  // Neutron visible E systs
426  else if( shift.first.find("NeutronEvis") != std::string::npos )
427  {
428  std::cerr << "Neutron uncert adding: -->" << shift.first << std::endl;
429  tempRet[kNeutronSys].push_back(shift.second);
430  continue;
431  }
432 
433  // Detector response systs
434  else if( shift.first.find("Cherenkov") != std::string::npos ||
435  shift.first.find("Light Level") != std::string::npos )
436  {
437  std::cerr << "Detector response adding: -->" << shift.first << std::endl;
438  tempRet[kBirksSys].push_back(shift.second);
439  continue;
440  }
441 
442  // Beam flux systs
443  else if( shift.first.find("ppfx") != std::string::npos )
444  {
445  std::cerr << "Beam flux adding: -->" << shift.first << std::endl;
446  tempRet[kBeamSys].push_back(shift.second);
447  continue;
448  }
449 
450  // Else Genie Xsec systs
451  else
452  {
453  std::cerr << "Neutrino cross sections adding: -->" << shift.first << std::endl;
454  tempRet[kGenieSys].push_back(shift.second);
455  continue;
456  }
457  }
458 
459  for (int i = 0; i < kNumSystTypes; ++i){
460  ret[systTypes[i]] = tempRet[i];
461  }
462 
463  return ret;
464 }
465 
467 {
468  TH1D* cosmicErrHist = cosmic->ToTH1(cosmic->Livetime(), kLivetime);
469  TH1D* cosmicShift = cosmic->ToTH1(cosmic->Livetime(), kLivetime);
470 
471  TFeldmanCousins fc(0.6827);
472  for(int i = 1; i<= cosmicErrHist->GetNbinsX(); i++){
473  double y = cosmicErrHist->GetBinContent(i);
474  double errup = 0;
475  double errlow = 0;
476  if(y!=0){
477  if ( y < 50 ) errlow = (y - fc.CalculateLowerLimit(y,0))/y;
478  else errlow = (sqrt(y))/y;
479  if ( y < 30 ) errup = (fc.CalculateUpperLimit(y,0)-y)/y;
480  else errup = (sqrt(y))/y;
481  }
482  if(sigma > 0)
483  cosmicShift->SetBinContent(i, cosmicShift->GetBinContent(i)*(1+sigma*errup));
484  if(sigma < 0)
485  cosmicShift->SetBinContent(i, cosmicShift->GetBinContent(i)*(1+sigma*errlow));
486  }
487 
488  Spectrum *ret = new Spectrum(cosmicShift, cosmic->POT(), cosmic->Livetime());
489  return ret;
490 }
T max(const caf::Proxy< T > &a, T b)
enum BeamMode kOrange
TSpline3 lo("lo", xlo, ylo, 12,"0")
enum BeamMode kRed
void drawBeamLabel(bool isFHC, double x=0.35)
set< int >::iterator it
std::map< std::string, std::vector< TH1 * > > SortSystHists(std::map< std::string, TH1 * > shifts)
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:148
T sqrt(T number)
Definition: d0nt_math.hpp:156
OStream cerr
Definition: OStream.cxx:7
Spectrum * ShiftedCosmics(Spectrum *cosmic, double sigma)
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:249
const Color_t kTotalMCErrorBandColor
Definition: Style.h:17
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
void documentHeader(bool IsNuMu)
const XML_Char * s
Definition: expat.h:262
bool IsNuMu(int pdgc)
Definition: PDGUtils.cxx:155
Current::Current_t curr
const int nbins
Definition: cellShifts.C:15
TSpline3 hi("hi", xhi, yhi, 18,"0")
void myBarChart(std::string xlabel, std::vector< std::string > ylabels, std::vector< std::vector< std::pair< double, double >>> errs, std::vector< std::pair< double, double >> sum, std::vector< std::string > legend, std::vector< int > colors)
int colors[6]
Definition: tools.h:1
Sum up livetimes from individual cosmic triggers.
std::vector< TH1 * > MakeRatios(TH1 *nom, std::vector< TH1 * > shifts)
Int_t col[ntarg]
Definition: Style.C:29
EExposureType
For use as an argument to Spectrum::ToTH1.
Definition: UtilsExt.h:35
void drawLabel(std::string s, double x=0.3, double y=0.93)
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
std::string systTypes[kNumSystTypes]
const double j
Definition: BetheBloch.cxx:29
std::vector< float > Spectrum
Definition: Constants.h:610
std::string name
TLatex * prelim
Definition: Xsec_final.C:133
double POT() const
Definition: Spectrum.h:227
static bool isFHC
double sigma(TH1F *hist, double percentile)
OStream cout
Definition: OStream.cxx:6
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
void drawSampleLabel(std::string s, double x=.915, double y=.9)
Sign::Sign_t sign
ifstream in
Definition: comparison.C:7
void MyPlotWithSystErrorBand(TH1 *&nom, std::vector< TH1 * > &ups, std::vector< TH1 * > &dns, int col, int errCol, float linewidth, bool newaxis)
enum BeamMode kViolet
TLatex * tex
Definition: f2_nu.C:499
TString FixSystName(TString mystring)
std::map< std::string, std::pair< double, double > > SumSysts(std::map< std::string, std::pair< double, double >> bars)
Flavors::Flavors_t flav
void preliminary()
const Color_t kTotalMCColor
Definition: Style.h:16
int systCols[kNumSystTypes]
T min(const caf::Proxy< T > &a, T b)
Double_t sum
Definition: plot.C:31
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
enum BeamMode kGreen
auto zero()
Definition: PMNS.cxx:47
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:230
double Integral(const Spectrum &s, double exp, double *err=0, EExposureType expotype=kPOT)
enum BeamMode string