NuSSystFunctions.h
Go to the documentation of this file.
1 /*
2  This file contains function systematic studies
3  In this file:
4 
5 
6  std::string StringFromInt(int i, bool withE20);
7  std::string StringFromDouble(double pot);
8 
9  void SqrtError(TH1* nom);
10  void ZeroError(TH1* nom);
11 
12  void AddErrorInQuadrature(TH1* nom, TH1* up, bool use_max);
13  void AddErrorInQuadrature(TH1* nom, TH1* up, TH1* down, bool use_max);
14 
15  TH1* GetNC(IDecomp* specs, double POT);
16  TH1* GetNC(IPrediction* specs, double POT);
17 
18  TH1* GetBG(IDecomp* specs, double POT);
19  TH1* GetBG(IPrediction* specs, double POT);
20 
21  void InitializeSystText(FILE* text, strings strs);
22 
23  void LoadMaps(TDirectory* dir,
24  std::map<std::string, IDecomp*> * nominal,
25  std::map<std::string, std::map<std::string, std::map<int, IDecomp*> > >* shifted);
26  void LoadMaps(TDirectory* dir,
27  std::map<std::string, PredictionNoExtrap*> * nominal,
28  std::map<std::string, std::map<std::string, std::map<int, PredictionNoExtrap*> > >* shifted);
29  void LoadMaps(TDirectory* dir,
30  std::map<std::string, PredictionSterile*> * nominal,
31  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > >* shifted);
32 
33  void PlotSyst(TH1* h0, TH1* hp1,
34  TDirectory* out, FILE* text, strings strs);
35  void PlotSyst(TH1* h0, TH1* hp1, TH1* hm1, TH1* hp2, TH1* hm2,
36  TDirectory* out, FILE* text, strings strs);
37 
38  void PlotSyst(TH1* nom, IDecomp* shifts,
39  TDirectory* out, FILE* text, strings strs,
40  double POT, bool NC = true, bool use_max);
41  void PlotSyst(TH1* nom, PredictionNoExtrap* shifts,
42  TDirectory* out, FILE* text, strings strs,
43  double POT, bool NC = true, bool use_max);
44  void PlotSyst(TH1* nom, PredictionSterile* shifts,
45  TDirectory* out, FILE* text, strings strs,
46  double POT, bool NC = true, bool use_max);
47 
48  void PlotSyst(TH1* nom, std::map<int, IDecomp*> shifts,
49  TDirectory* out, FILE* text, strings strs,
50  double POT, bool NC = true, bool use_max);
51  void PlotSyst(TH1* nom, std::map<int, PredictionNoExtrap*> shifts,
52  TDirectory* out, FILE* text, strings strs,
53  double POT, bool NC = true, bool use_max);
54  void PlotSyst(TH1* nom, std::map<int, PredictionSterile*> shifts,
55  TDirectory* out, FILE* text, strings strs,
56  double POT, bool NC = true, bool use_max);
57 
58  void PlotSystBand(TH1* h, TDirectory* out, FILE* text, strings strs);
59 
60  void SaveMaps(
61  TDirectory* out,
62  std::map<std::string, IDecomp*> decomps_nominal,
63  std::map<std::string, std::map<std::string, std::map<int, IDecomp*> > > decomps_shifted,
64  std::map<std::string, PredictionNoExtrap*> predsNE_nominal,
65  std::map<std::string, std::map<std::string, std::map<int, PredictionNoExtrap*> > > predsNE_shifted,
66  std::map<std::string, PredictionSterile*> predsSt_nominal,
67  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > > predsSt_shifted
68  );
69 */
70 
71 #pragma once
72 
73 #include <map>
74 #include <iostream>
75 #include <string>
76 #include <vector>
77 
78 #include "CAFAna/Analysis/Calcs.h"
79 #include "CAFAna/Analysis/Style.h"
81 #include "CAFAna/Decomp/IDecomp.h"
87 
88 #include "OscLib/OscCalcPMNSOpt.h"
89 #include "OscLib/OscCalcSterile.h"
90 
91 #include "TCanvas.h"
92 #include "TDirectory.h"
93 #include "TH1.h"
94 #include "TKey.h"
95 #include "TLatex.h"
96 #include "TLegend.h"
97 
98 namespace ana
99 {
100  /// A helper structure to contain a group of string for plotting
101  struct strings
102  {
111  };
112 
113  //----------------------------------------------------------------------------
114  /// Helper function that takes an integer and turns it into a string
115  std::string StringFromInt(int i, bool withE20)
116  {
117  char pot[3];
118  sprintf(pot, "%d", i);
119  std::string potStr = std::string(pot);
120  if(withE20) { potStr += " x 10^{20} POT"; }
121 
122  return potStr;
123  }
124 
125  //----------------------------------------------------------------------------
126  /// Helper function that takes a double POT value and turns it into a string
127  /// The output has the format "# x 10^{20} POT"
128  std::string StringFromDouble(double pot)
129  {
130  double dropE20 = pot/1e20;
131  int i = (int)dropE20;
132  if(dropE20 > (double)i + 0.5) { ++i; }
133 
134  return StringFromInt(i, true);
135  }
136 
137  //----------------------------------------------------------------------------
138  /// Set the errors of the input histogram to be their square roots
139  void SqrtError(TH1* nom)
140  {
141  for(int i = 0, n = nom->GetNbinsX(); i <= n; ++i) {
142  nom->SetBinError(i, sqrt(nom->GetBinError(i)));
143  }
144  return;
145  }
146 
147  //----------------------------------------------------------------------------
148  /// Set the errors of the input histogram to 0
149  void ZeroError(TH1* nom)
150  {
151  for(int i = 0, n = nom->GetNbinsX(); i <= n; ++i) {
152  nom->SetBinError(i, 0.);
153  }
154  return;
155  }
156 
157  //----------------------------------------------------------------------------
158  /// Combine errors in quadrature, setting error bins in nominal spectra
159  /// Leaves errors squared -- square root must be applied later
160  /// \param nom The nominal histogram
161  /// \param up The shifted spectra
162  /// \param use_max Only keep the max error, don't combine in quadrature
163  void AddErrorInQuadrature(TH1* nom, TH1* up, bool use_max)
164  {
165  double error;
166  for(int i = 0, n = nom->GetNbinsX(); i <= n; ++i) {
167  error = std::abs(up->GetBinContent(i) - nom->GetBinContent(i));
168  double error_sq = error*error;
169  double err_prev = std::abs(nom->GetBinError(i));
170  if(use_max) {
171  nom->SetBinError(i, std::max(err_prev, error_sq));
172  }
173  else {
174  nom->SetBinError(i, err_prev + error_sq);
175  }
176  }
177 
178  return;
179  }
180 
181  //----------------------------------------------------------------------------
182  /// Combine errors in quadrature, setting error bins in nominal spectra
183  /// Takes the larger difference from nominal from either up or down shift
184  /// Leaves errors squared -- square root must be applied later
185  /// \param nom The nominal histogram
186  /// \param up The +1 sigma shifted spectra
187  /// \param down The -1 sigma shifted spectra
188  /// \param use_max Only keep the max error, don't combine in quadrature
189  void AddErrorInQuadrature(TH1* nom, TH1* up, TH1* down, bool use_max)
190  {
191  double error;
192  for(int i = 0, n = nom->GetNbinsX(); i <= n; ++i) {
193  error = std::max(std::abs(up ->GetBinContent(i) - nom ->GetBinContent(i)),
194  std::abs(nom->GetBinContent(i) - down->GetBinContent(i)));
195  double error_sq = error*error;
196  double err_prev = std::abs(nom->GetBinError(i));
197  if(use_max) {
198  nom->SetBinError(i, std::max(err_prev, error_sq));
199  }
200  else {
201  nom->SetBinError(i, err_prev + error_sq);
202  }
203  }
204 
205  return;
206  }
207 
208  //----------------------------------------------------------------------------
209  /// Get the NC spectrum from a decomposition at a specific POT value
210  /// Sets errors to 0
211  TH1* GetNC(IDecomp* specs, double POT)
212  {
213  Spectrum s = specs->NCTotalComponent();
214  TH1* ret = s.ToTH1(POT);
215  ZeroError(ret);
216  return ret;
217  }
218 
219  //----------------------------------------------------------------------------
220  /// Get the NC spectrum from a prediction object at a specific POT value
221  /// Sets errors to 0
222  TH1* GetNC(IPrediction* specs, double POT)
223  {
225  Spectrum s = specs->PredictComponent(
227  );
228  TH1* ret = s.ToTH1(POT);
229  ZeroError(ret);
230  return ret;
231  }
232 
233  //----------------------------------------------------------------------------
234  /// Get the background spectrum from a decomposition at a specific POT value
235  /// Sets errors to 0
236  TH1* GetBG(IDecomp* specs, double POT)
237  {
238  Spectrum s = specs->NueComponent();
239  s += specs->AntiNueComponent();
240  s += specs->NumuComponent();
241  s += specs->AntiNumuComponent();
242  TH1* ret = s.ToTH1(POT);
243  ZeroError(ret);
244  return ret;
245  }
246 
247  //----------------------------------------------------------------------------
248  /// Get the background spectrum from a prediction object at a specific POT value
249  /// Sets errors to 0
250  TH1* GetBG(IPrediction* specs, double POT)
251  {
253  Spectrum s = specs->PredictComponent(
255  );
256  TH1* ret = s.ToTH1(POT);
257  ZeroError(ret);
258  return ret;
259  }
260 
261  //----------------------------------------------------------------------------
262  /// Print some initial text about a systematic -- the systematic type and sample
263  void InitializeSystText(FILE* text, strings strs)
264  {
265  fprintf(text, "%s %s\n\n", strs.fSystType.c_str(), strs.fSample.c_str());
266 
267  return;
268  }
269 
270  //----------------------------------------------------------------------------
271  /// Load nominal and shifted decompositions from a TDirectory
272  /// The maps to contain these objects must already exist
273  void LoadMaps(TDirectory* dir,
274  std::map<std::string, IDecomp*> * nominal,
275  std::map<std::string, std::map<std::string, std::map<int, IDecomp*> > >* shifted)
276  {
277  TDirectory* tmp = gDirectory;
278  TDirectory* keyDir = gDirectory;
279 
280  // Loop over all keys in this directory
281  dir->cd();
282  TIter nextkey(dir->GetListOfKeys());
283  TKey* key;
284  TKey* oldkey = 0;
285 
286  std::string sep = "__"; // This is the token separator
287  std::string cafanaType = "decomp";
288  std::string sampleLabel = "";
289  std::string systLabel = "";
290  int sigma = 0;
291 
292  while((key = (TKey*)nextkey())) {
293  // Keep only the highest cycle number for each key
294  if(oldkey && !strcmp(oldkey->GetName(), key->GetName())) { continue; }
295  if(key->ReadObj()->IsFolder()) {
296  keyDir = (TDirectory*)key->ReadObj();
297  }
298  else {
299  continue;
300  }
301 
302  std::string name(key->GetName());
303  std::vector<std::string> tokens;
304  std::size_t pos = 0, prev = 0;
305  while(pos != std::string::npos) { // Loop over the TObject name
306  pos = name.find(sep, prev); // Find the next separator
307 
308  if(pos != std::string::npos) {
309  tokens.push_back(name.substr(prev, pos-prev));
310  prev = pos+2; // Move past the current separator
311  }
312  else {
313  tokens.push_back(name.substr(prev, name.length()-prev));
314  }
315  }
316 
317  // Make sure the object we found is the type we want
318  if(cafanaType.compare(tokens[0]) != 0) { continue; }
319 
320  sampleLabel = tokens[1];
321  if(tokens.size() > 2) {
322  systLabel = tokens[2];
323  sigma = std::stoi(tokens[3]);
324 
325  (*shifted)[sampleLabel][systLabel][sigma] = LoadFrom<IDecomp>(keyDir).release();
326  }
327  else {
328  (*nominal)[sampleLabel] = LoadFrom<IDecomp>(keyDir).release();
329  }
330  }
331 
332  tmp->cd();
333  return;
334  }
335 
336  //----------------------------------------------------------------------------
337  /// Load nominal and shifted PredictionNoExtrap objects from a TDirectory
338  /// The maps to contain these objects must already exist
339  void LoadMaps(TDirectory* dir,
340  std::map<std::string, PredictionNoExtrap*> * nominal,
341  std::map<std::string, std::map<std::string, std::map<int, PredictionNoExtrap*> > >* shifted)
342  {
343  TDirectory* tmp = gDirectory;
344  TDirectory* keyDir = gDirectory;
345 
346  dir->cd();
347  TIter nextkey(dir->GetListOfKeys());
348  TKey* key;
349  TKey* oldkey = 0;
350 
351  std::string sep = "__";
352  std::string cafanaType = "predNE";
353  std::string sampleLabel = "";
354  std::string systLabel = "";
355  int sigma = 0;
356 
357  while((key = (TKey*)nextkey())) {
358  if(oldkey && !strcmp(oldkey->GetName(), key->GetName())) { continue; }
359  if(key->ReadObj()->IsFolder()) {
360  keyDir = (TDirectory*)key->ReadObj();
361  }
362  else {
363  continue;
364  }
365 
366  std::string name(key->GetName());
367  std::vector<std::string> tokens;
368  std::size_t pos = 0, prev = 0;
369  while(pos != std::string::npos) {
370  pos = name.find(sep, prev);
371 
372  if(pos != std::string::npos) {
373  tokens.push_back(name.substr(prev, pos-prev));
374  prev = pos+2;
375  }
376  else {
377  tokens.push_back(name.substr(prev, name.length()-prev));
378  }
379  }
380 
381  if(cafanaType.compare(tokens[0]) != 0) { continue; }
382 
383  sampleLabel = tokens[1];
384  if(tokens.size() > 2) {
385  systLabel = tokens[2];
386  sigma = std::stoi(tokens[3]);
387 
388  (*shifted)[sampleLabel][systLabel][sigma] = (PredictionNoExtrap*)PredictionExtrap::LoadFrom(keyDir).release();
389  }
390  else {
391  (*nominal)[sampleLabel] = (PredictionNoExtrap*)PredictionExtrap::LoadFrom(keyDir).release();
392  }
393  }
394 
395  tmp->cd();
396  return;
397  }
398 
399  //----------------------------------------------------------------------------
400  /// Load nominal and shifted PredictionSterile objects from a TDirectory
401  /// The maps to contain these objects must already exist
402  void LoadMaps(TDirectory* dir,
403  std::map<std::string, PredictionSterile*> * nominal,
404  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > >* shifted)
405  {
406  TDirectory* tmp = gDirectory;
407  TDirectory* keyDir = gDirectory;
408 
409  dir->cd();
410  TIter nextkey(dir->GetListOfKeys());
411  TKey* key;
412  TKey* oldkey = 0;
413 
414  std::string sep = "__";
415  std::string cafanaType = "predSt";
416  std::string sampleLabel = "";
417  std::string systLabel = "";
418  int sigma = 0;
419 
420  while((key = (TKey*)nextkey())) {
421  if(oldkey && !strcmp(oldkey->GetName(), key->GetName())) { continue; }
422  if(key->ReadObj()->IsFolder()) {
423  keyDir = (TDirectory*)key->ReadObj();
424  }
425  else {
426  continue;
427  }
428 
429  std::string name(key->GetName());
430  std::vector<std::string> tokens;
431  std::size_t pos = 0, prev = 0;
432  while(pos != std::string::npos) {
433  pos = name.find(sep, prev);
434 
435  if(pos != std::string::npos) {
436  tokens.push_back(name.substr(prev, pos-prev));
437  prev = pos+2;
438  }
439  else {
440  tokens.push_back(name.substr(prev, name.length()-prev));
441  }
442  }
443 
444  if(cafanaType.compare(tokens[0]) != 0) { continue; }
445 
446  sampleLabel = tokens[1];
447  if(tokens.size() > 2) {
448  systLabel = tokens[2];
449  sigma = std::stoi(tokens[3]);
450 
451  (*shifted)[sampleLabel][systLabel][sigma] = PredictionSterile::LoadFrom(keyDir).release();
452  }
453  else {
454  (*nominal)[sampleLabel] = PredictionSterile::LoadFrom(keyDir).release();
455  }
456  }
457 
458  tmp->cd();
459  return;
460  }
461 
462  //----------------------------------------------------------------------------
463  /// Plot a nominal spectrum and a shifted spectra
464  /// This is for an on/off sytle systematic
465  /// \param h The nominal spectrum (systematic off)
466  /// \param hp1 The shifted spectrum (systematic on)
467  /// \param out The directory/file for root output
468  /// \param text The file for text output
469  /// \param strs Contains strings for pretty plots
470  void PlotSyst(TH1* h, TH1* hp1, TDirectory* out, FILE* text, strings strs)
471  {
472  // Construct histogram names and titles
473  std::string xlabel = strs.fXLabel;
474  std::string ylabel = "Events / " + strs.fPOT + " / 0.25 GeV";
475  std::string title = strs.fDet + " " + strs.fComponent +
476  " Error for " + strs.fSystL + " Systematic";
477  std::string fullTitle = ";" + xlabel + ";" + ylabel;
478 // std::string fullTitle = title + ";" + xlabel + ";" + ylabel;
479  std::string name = "c" + strs.fComponent + strs.fDet +
480  strs.fSystS + strs.fSample;
481  std::string ylabelRat = "Shifted / Nominal";
482  std::string fullTitleRat = ";" + xlabel + ";" + ylabelRat;
483 
484  TH1* hRat = (TH1*)h->Clone();
485  TH1* hp1Rat = (TH1*)hp1->Clone();
486 
487  int nBins = h->GetNbinsX();
488  double percent_diff = 100. * std::abs(hp1->Integral(1, nBins) - h->Integral(1, nBins))/
489  h->Integral(1, nBins);
490  fprintf(text, "%s %s %.2s: +/-%6.4f%%\n",
491  strs.fSystS.c_str(), strs.fComponent.c_str(), strs.fDet.c_str(),
492  percent_diff);
493 
494  for(int i = 0; i <= nBins+1; ++i) {
495  hRat->SetBinContent(i, 1.);
496  }
497  hp1Rat->Divide(h);
498 
499  double maxval = h->GetBinContent(h->GetMaximumBin());
500  maxval = std::max(maxval, hp1->GetBinContent(hp1->GetMaximumBin()));
501 
502  double minval = h->GetBinContent(h->GetMinimumBin());
503  minval = std::min(minval, hp1->GetBinContent(hp1->GetMinimumBin()));
504 
505  if(maxval < minval) { std::swap(maxval, minval); }
506 
507  double maxvalRat = hRat->GetBinContent(hRat->GetMaximumBin());
508  maxvalRat = std::max(maxvalRat, hp1Rat->GetBinContent(hp1Rat->GetMaximumBin()));
509 
510  double minvalRat = hRat->GetBinContent(hRat->GetMinimumBin());
511  for(int i = 1; i <= nBins; ++i) {
512  if(hp1Rat->GetBinContent(i) > 0.) {
513  minvalRat = std::min(minvalRat, hp1Rat->GetBinContent(i));
514  }
515  }
516 
517  if(std::abs(maxvalRat - 1.) > std::abs(1. - minvalRat)) {
518  minvalRat = 1 - std::abs(maxvalRat - 1.);
519  maxvalRat = 1 + std::abs(maxvalRat - 1.);
520  }
521  else {
522  minvalRat = 1 - std::abs(1. - minvalRat);
523  maxvalRat = 1 + std::abs(1. - minvalRat);
524  }
525 
526  CenterTitles(h);
527  h->SetLineColor(kBlack);
528  h->SetMaximum(1.1*maxval);
529  h->SetMinimum(0);
530  h->SetTitle(fullTitle.c_str());
531 
532  CenterTitles(hp1);
533  hp1->SetLineColor(kRed);
534  hp1->SetLineStyle(2);
535  hp1->SetMaximum(1.1*maxval);
536  hp1->SetMinimum(0);
537  hp1->SetTitle(fullTitle.c_str());
538 
539  CenterTitles(hRat);
540  hRat->SetLineColor(kBlack);
541  hRat->SetMaximum(1.1*maxvalRat);
542  hRat->SetMinimum(0.9*minvalRat);
543  hRat->SetTitle(fullTitleRat.c_str());
544 
545  CenterTitles(hp1Rat);
546  hp1Rat->SetLineColor(kRed);
547  hp1Rat->SetLineStyle(2);
548  hp1Rat->SetMaximum(1.1*maxvalRat);
549  hp1Rat->SetMinimum(0.9*minvalRat);
550  hp1Rat->SetTitle(fullTitleRat.c_str());
551 
552  double xL = 0.6, xR = 0.85;
553  double yB = 0.6, yT = 0.85;
554  TLegend* leg = new TLegend(xL, yB, xR, yT);
555  leg->AddEntry(h, "Nominal", "l");
556  leg->AddEntry(hp1, "Shifted", "l");
557 
558  TCanvas* c = new TCanvas(name.c_str(), title.c_str(), 800, 800);
559  TPad* pSpecs = new TPad("pSpecs", "", 0., 0.375, 1., 1.);
560  TPad* pRatio = new TPad("pRatio", "", 0., 0., 1., 0.375);
561  pSpecs->Draw();
562  pRatio->Draw();
563 
564  pSpecs->cd();
565  h->Draw("hist");
566  hp1->Draw("hist same");
567  leg->Draw();
568  Simulation();
569 
570  pRatio->cd();
571  hRat->Draw("hist");
572  hp1Rat->Draw("hist same");
573 
574  c->Update();
575  out->WriteTObject(c);
576  return;
577  }
578 
579  //----------------------------------------------------------------------------
580  /// Plot a nominal spectrum and a shifted spectra
581  /// This is for a systematic with multiple sigma levels
582  /// \param h The nominal spectrum
583  /// \param hp1 The shifted spectrum, +1 sigma
584  /// \param hm1 The shifted spectrum, -1 sigma
585  /// \param hp2 The shifted spectrum, +2 sigma
586  /// \param hm2 The shifted spectrum, -2 sigma
587  /// \param out The directory/file for root output
588  /// \param text The file for text output
589  /// \param strs Contains strings for pretty plots
590  void PlotSyst(TH1* h, TH1* hp1, TH1* hm1, TH1* hp2, TH1* hm2,
591  TDirectory* out, FILE* text, strings strs)
592  {
593  std::string xlabel = strs.fXLabel;
594  std::string ylabel = "Events / " + strs.fPOT + " / 0.25 GeV";
595  std::string title = strs.fDet + " " + strs.fComponent +
596  " Error for " + strs.fSystL + " Systematic";
597  std::string fullTitle = ";" + xlabel + ";" + ylabel;
598 // std::string fullTitle = title + ";" + xlabel + ";" + ylabel;
599  std::string name = "c" + strs.fComponent + strs.fDet +
600  strs.fSystS + strs.fSample;
601  std::string ylabelRat = "Shifted / Nominal";
602  std::string fullTitleRat = ";" + xlabel + ";" + ylabelRat;
603 
604  TH1* hRat = (TH1*)h->Clone();
605  TH1* hp1Rat = (TH1*)hp1->Clone();
606  TH1* hm1Rat = (TH1*)hm1->Clone();
607  TH1* hp2Rat = (TH1*)hp2->Clone();
608  TH1* hm2Rat = (TH1*)hm2->Clone();
609 
610  int nBins = h->GetNbinsX();
611  double nom_evts = h->Integral(1, nBins);
612  double evts_diff = 0.;
613  double percent_diff = 100.;
614  double percent_diff_p = 100.;
615  double percent_diff_m = 100.;
616  for(int i = 0; i <= nBins+1; ++i) {
617  hRat->SetBinContent(i, 1.);
618 
619  if(i != 0 && i != nBins+1) {
620  evts_diff += std::max(std::abs(hp1->GetBinContent(i) - h->GetBinContent(i)),
621  std::abs(hm1->GetBinContent(i) - h->GetBinContent(i)));
622  }
623  }
624  percent_diff *= (evts_diff/nom_evts);
625  percent_diff_p *= ((hp1->Integral(1, nBins) - nom_evts)/nom_evts);
626  percent_diff_m *= ((nom_evts - hm1->Integral(1, nBins))/nom_evts);
627 
628  fprintf(text, "%s %.2s %.2s: +/-%6.4f%%, +%6.4f%%, -%6.4f%%\n",
629  strs.fSystS.c_str(), strs.fComponent.c_str(), strs.fDet.c_str(),
630  percent_diff, percent_diff_p, percent_diff_m);
631 
632  hp1Rat->Divide(h);
633  hm1Rat->Divide(h);
634  hp2Rat->Divide(h);
635  hm2Rat->Divide(h);
636 
637  double maxval = h->GetBinContent(h->GetMaximumBin());
638  maxval = std::max(maxval, hp1->GetBinContent(hp1->GetMaximumBin()));
639  maxval = std::max(maxval, hm1->GetBinContent(hm1->GetMaximumBin()));
640  maxval = std::max(maxval, hp2->GetBinContent(hp2->GetMaximumBin()));
641  maxval = std::max(maxval, hm2->GetBinContent(hm2->GetMaximumBin()));
642 
643  double minval = h->GetBinContent(h->GetMinimumBin());
644  minval = std::min(minval, hp1->GetBinContent(hp1->GetMinimumBin()));
645  minval = std::min(minval, hm1->GetBinContent(hm1->GetMinimumBin()));
646  minval = std::min(minval, hp2->GetBinContent(hp2->GetMinimumBin()));
647  minval = std::min(minval, hm2->GetBinContent(hm2->GetMinimumBin()));
648 
649  if(maxval < minval) { std::swap(maxval, minval); }
650 
651  double maxvalRat = hRat->GetBinContent(hRat->GetMaximumBin());
652  maxvalRat = std::max(maxvalRat, hp1Rat->GetBinContent(hp1Rat->GetMaximumBin()));
653  maxvalRat = std::max(maxvalRat, hm1Rat->GetBinContent(hm1Rat->GetMaximumBin()));
654  maxvalRat = std::max(maxvalRat, hp2Rat->GetBinContent(hp2Rat->GetMaximumBin()));
655  maxvalRat = std::max(maxvalRat, hm2Rat->GetBinContent(hm2Rat->GetMaximumBin()));
656 
657  double minvalRat = hRat->GetBinContent(hRat->GetMinimumBin());
658  for(int i = 1; i <= nBins; ++i) {
659  if(hp1Rat->GetBinContent(i) > 0.) {
660  minvalRat = std::min(minvalRat, hp1Rat->GetBinContent(i));
661  }
662  if(hm1Rat->GetBinContent(i) > 0.) {
663  minvalRat = std::min(minvalRat, hm1Rat->GetBinContent(i));
664  }
665  if(hp2Rat->GetBinContent(i) > 0.) {
666  minvalRat = std::min(minvalRat, hp2Rat->GetBinContent(i));
667  }
668  if(hm2Rat->GetBinContent(i) > 0.) {
669  minvalRat = std::min(minvalRat, hm2Rat->GetBinContent(i));
670  }
671  }
672 
673  if(std::abs(maxvalRat - 1.) > std::abs(1. - minvalRat)) {
674  minvalRat = 1 - std::abs(maxvalRat - 1.);
675  maxvalRat = 1 + std::abs(maxvalRat - 1.);
676  }
677  else {
678  minvalRat = 1 - std::abs(1. - minvalRat);
679  maxvalRat = 1 + std::abs(1. - minvalRat);
680  }
681 
682  CenterTitles(h);
683  h->SetLineColor(kBlack);
684  h->SetMaximum(1.1*maxval);
685  h->SetMinimum(0);
686  h->SetTitle(fullTitle.c_str());
687 
688  CenterTitles(hp1);
689  hp1->SetLineColor(kRed);
690  hp1->SetLineStyle(2);
691  hp1->SetMaximum(1.1*maxval);
692  hp1->SetMinimum(0);
693  hp1->SetTitle(fullTitle.c_str());
694 
695  CenterTitles(hm1);
696  hm1->SetLineColor(kRed);
697  hm1->SetLineStyle(2);
698  hm1->SetMaximum(1.1*maxval);
699  hm1->SetMinimum(0);
700  hm1->SetTitle(fullTitle.c_str());
701 
702  CenterTitles(hp2);
703  hp2->SetLineColor(kBlue);
704  hp2->SetLineStyle(2);
705  hp2->SetMaximum(1.1*maxval);
706  hp2->SetMinimum(0);
707  hp2->SetTitle(fullTitle.c_str());
708 
709  CenterTitles(hm2);
710  hm2->SetLineColor(kBlue);
711  hm2->SetLineStyle(2);
712  hm2->SetMaximum(1.1*maxval);
713  hm2->SetMinimum(0);
714  hm2->SetTitle(fullTitle.c_str());
715 
716  CenterTitles(hRat);
717  hRat->SetLineColor(kBlack);
718  hRat->SetMaximum(1.1*maxvalRat);
719  hRat->SetMinimum(0.9*minvalRat);
720  hRat->SetTitle(fullTitleRat.c_str());
721 
722  CenterTitles(hp1Rat);
723  hp1Rat->SetLineColor(kRed);
724  hp1Rat->SetLineStyle(2);
725  hp1Rat->SetMaximum(1.1*maxvalRat);
726  hp1Rat->SetMinimum(0.9*minvalRat);
727  hp1Rat->SetTitle(fullTitleRat.c_str());
728 
729  CenterTitles(hm1Rat);
730  hm1Rat->SetLineColor(kRed);
731  hm1Rat->SetLineStyle(2);
732  hm1Rat->SetMaximum(1.1*maxvalRat);
733  hm1Rat->SetMinimum(0.9*minvalRat);
734  hm1Rat->SetTitle(fullTitleRat.c_str());
735 
736  CenterTitles(hp2Rat);
737  hp2Rat->SetLineColor(kBlue);
738  hp2Rat->SetLineStyle(2);
739  hp2Rat->SetMaximum(1.1*maxvalRat);
740  hp2Rat->SetMinimum(0.9*minvalRat);
741  hp2Rat->SetTitle(fullTitleRat.c_str());
742 
743  CenterTitles(hm2Rat);
744  hm2Rat->SetLineColor(kBlue);
745  hm2Rat->SetLineStyle(2);
746  hm2Rat->SetMaximum(1.1*maxvalRat);
747  hm2Rat->SetMinimum(0.9*minvalRat);
748  hm2Rat->SetTitle(fullTitleRat.c_str());
749 
750  double xL = 0.6, xR = 0.85;
751  double yB = 0.6, yT = 0.85;
752  TLegend* leg = new TLegend(xL, yB, xR, yT);
753  leg->AddEntry(h, "Nominal", "l");
754  leg->AddEntry(hp1, "+/- 1 #sigma", "l");
755  leg->AddEntry(hp2, "+/- 2 #sigma", "l");
756 
757  TCanvas* c = new TCanvas(name.c_str(), title.c_str(), 800, 800);
758  TPad* pSpecs = new TPad("pSpecs", "", 0., 0.375, 1., 1.);
759  TPad* pRatio = new TPad("pRatio", "", 0., 0., 1., 0.375);
760  pSpecs->Draw();
761  pRatio->Draw();
762 
763  pSpecs->cd();
764  h->Draw("hist");
765  hp2->Draw("hist same");
766  hm2->Draw("hist same");
767  hp1->Draw("hist same");
768  hm1->Draw("hist same");
769  leg->Draw();
770  Simulation();
771 
772  pRatio->cd();
773  hRat->Draw("hist");
774  hp2Rat->Draw("hist same");
775  hm2Rat->Draw("hist same");
776  hp1Rat->Draw("hist same");
777  hm1Rat->Draw("hist same");
778 
779  c->Update();
780  out->WriteTObject(c);
781  return;
782  }
783 
784  //----------------------------------------------------------------------------
785  /// Plot ND syst shifted and nominal spectra
786  /// This is for plotting an on/off style shift
787  void PlotSyst(TH1* nom, IDecomp* shifts,
788  TDirectory* out, FILE* text, strings strs,
789  double POT, bool NC, bool use_max)
790  {
791  TH1* hp1 = (NC ? GetNC(shifts, POT) : GetBG(shifts, POT));
792 
793  AddErrorInQuadrature(nom, hp1, use_max);
794  PlotSyst(nom, hp1, out, text, strs);
795  return;
796  }
797 
798  //----------------------------------------------------------------------------
799  /// Plot FD syst shifted and nominal spectra
800  /// Should be able to change the input to IPrediction
801  /// This is for plotting an on/off style shift
802  void PlotSyst(TH1* nom, PredictionNoExtrap* shifts,
803  TDirectory* out, FILE* text, strings strs,
804  double POT, bool NC, bool use_max)
805  {
806  TH1* hp1 = (NC ? GetNC(shifts, POT) : GetBG(shifts, POT));
807 
808  AddErrorInQuadrature(nom, hp1, use_max);
809  PlotSyst(nom, hp1, out, text, strs);
810  return;
811  }
812 
813  //----------------------------------------------------------------------------
814  /// Plot FD syst shifted and nominal spectra
815  /// Should be able to change the input to IPrediction and make this redundant
816  /// This is for plotting an on/off style shift
817  void PlotSyst(TH1* nom, PredictionSterile* shifts,
818  TDirectory* out, FILE* text, strings strs,
819  double POT, bool NC, bool use_max)
820  {
821  TH1* hp1 = (NC ? GetNC(shifts, POT) : GetBG(shifts, POT));
822 
823  AddErrorInQuadrature(nom, hp1, use_max);
824  PlotSyst(nom, hp1, out, text, strs);
825  return;
826  }
827 
828  //----------------------------------------------------------------------------
829  /// Plot ND syst shifted and nominal spectra
830  /// This is for plotting a syst with multiple sigma levels
831  void PlotSyst(TH1* nom, std::map<int, IDecomp*> shifts,
832  TDirectory* out, FILE* text, strings strs,
833  double POT, bool NC, bool use_max)
834  {
835  TH1* hp1 = (NC ? GetNC(shifts[+1], POT) : GetBG(shifts[+1], POT));
836  TH1* hm1 = (NC ? GetNC(shifts[-1], POT) : GetBG(shifts[-1], POT));
837  TH1* hp2 = (NC ? GetNC(shifts[+2], POT) : GetBG(shifts[+2], POT));
838  TH1* hm2 = (NC ? GetNC(shifts[-2], POT) : GetBG(shifts[-2], POT));
839 
840  AddErrorInQuadrature(nom, hp1, hm1, use_max);
841  PlotSyst(nom, hp1, hm1, hp2, hm2, out, text, strs);
842  return;
843  }
844 
845  //----------------------------------------------------------------------------
846  /// Plot FD syst shifted and nominal spectra
847  /// Should be able to change the input to IPrediction
848  /// This is for plotting a syst with multiple sigma levels
849  void PlotSyst(TH1* nom, std::map<int, PredictionNoExtrap*> shifts,
850  TDirectory* out, FILE* text, strings strs,
851  double POT, bool NC, bool use_max)
852  {
853  TH1* hp1 = (NC ? GetNC(shifts[+1], POT) : GetBG(shifts[+1], POT));
854  TH1* hm1 = (NC ? GetNC(shifts[-1], POT) : GetBG(shifts[-1], POT));
855  TH1* hp2 = (NC ? GetNC(shifts[+2], POT) : GetBG(shifts[+2], POT));
856  TH1* hm2 = (NC ? GetNC(shifts[-2], POT) : GetBG(shifts[-2], POT));
857 
858  AddErrorInQuadrature(nom, hp1, hm1, use_max);
859  PlotSyst(nom, hp1, hm1, hp2, hm2, out, text, strs);
860  return;
861  }
862 
863  //----------------------------------------------------------------------------
864  /// Plot FD syst shifted and nominal spectra
865  /// Should be able to change the input to IPrediction and make this redundant
866  /// This is for plotting a syst with multiple sigma levels
867  void PlotSyst(TH1* nom, std::map<int, PredictionSterile*> shifts,
868  TDirectory* out, FILE* text, strings strs,
869  double POT, bool NC, bool use_max)
870  {
871  TH1* hp1 = (NC ? GetNC(shifts[+1], POT) : GetBG(shifts[+1], POT));
872  TH1* hm1 = (NC ? GetNC(shifts[-1], POT) : GetBG(shifts[-1], POT));
873  TH1* hp2 = (NC ? GetNC(shifts[+2], POT) : GetBG(shifts[+2], POT));
874  TH1* hm2 = (NC ? GetNC(shifts[-2], POT) : GetBG(shifts[-2], POT));
875 
876  AddErrorInQuadrature(nom, hp1, hm1, use_max);
877  PlotSyst(nom, hp1, hm1, hp2, hm2, out, text, strs);
878  return;
879  }
880 
881  //----------------------------------------------------------------------------
882  /// Plot the systematic error band around a nominal spectrum
883  /// The errors of the input histogram should be squared
884  void PlotSystBand(TH1* h, TDirectory* out, FILE* text, strings strs)
885  {
886  // Fix the errors by taking the square root!
887  SqrtError(h);
888 
889  std::string xlabel = strs.fXLabel;
890  std::string ylabel = "Events / " + strs.fPOT + " / 0.25 GeV";
891  std::string title = strs.fDet + " " + strs.fComponent +
892  " Error Band for All " + strs.fSystType + " Systematics";
893  std::string fullTitle = ";" + xlabel + ";" + ylabel;
894 // std::string fullTitle = title + ";" + xlabel + ";" + ylabel;
895  std::string name = "c" + strs.fComponent + strs.fDet +
896  strs.fSystType + "Systs" + strs.fSample;
897  std::string ylabelRat = "Shifted / Nominal";
898  std::string fullTitleRat = ";" + xlabel + ";" + ylabelRat;
899 
900  int nBins = h->GetNbinsX();
901  double maxval = 0.;
902  double nom_evts = h->Integral(1, nBins);
903  double evts_diff = 0.;
904  double percent_diff = 100.;
905  for(int i = 1; i <= nBins; ++i) {
906  maxval = std::max(maxval, h->GetBinContent(i) + h->GetBinError(i));
907  evts_diff += std::abs(h->GetBinError(i));
908  }
909  percent_diff *= (evts_diff/nom_evts);
910 
911  fprintf(text, "\nAll systs %.2s %.2s: +/-%6.4f%%\n",
912  strs.fComponent.c_str(), strs.fDet.c_str(),
913  percent_diff);
914 
915  CenterTitles(h);
916  h->SetLineColor(kTotalMCColor);
917  h->SetMaximum(1.1*maxval);
918  h->SetMinimum(0);
919  h->SetTitle(fullTitle.c_str());
920 
921  TH1* hErr = (TH1*)h->Clone();
922  CenterTitles(hErr);
923  hErr->SetFillColor(kTotalMCColor-9);
924  hErr->SetMaximum(1.1*maxval);
925  hErr->SetMinimum(0);
926 
927  double maxvalRat = 1.;
928  TH1* hRat = (TH1*)h->Clone();
929  TH1* hUpp = (TH1*)h->Clone();
930  TH1* hDwn = (TH1*)h->Clone();
931  for(int i = 1, n = hRat->GetNbinsX(); i <= n; ++i) {
932  double binVal = hRat->GetBinContent(i);
933  double errVal = hRat->GetBinError(i);
934  if(binVal <= 0.) {
935  binVal = 1.;
936  errVal = 0.;
937  }
938 
939  hRat->SetBinContent(i, 1);
940  hUpp->SetBinContent(i, (binVal + errVal)/binVal);
941  hDwn->SetBinContent(i, (binVal - errVal)/binVal);
942 
943  double diff = (binVal + errVal)/binVal;
944  if(diff > maxvalRat) {
945  maxvalRat = diff;
946  }
947  }
948  double minvalRat = 1. - (maxvalRat-1.);
949 
950  CenterTitles(hRat);
951  hRat->SetLineColor(kBlack);
952  hRat->SetMaximum(1.05*maxvalRat);
953  hRat->SetMinimum(0.95*minvalRat);
954  hRat->SetTitle(fullTitleRat.c_str());
955 
956  CenterTitles(hUpp);
957  hUpp->SetLineColor(kRed);
958  hUpp->SetLineStyle(2);
959  hUpp->SetMaximum(1.05*maxvalRat);
960  hUpp->SetMinimum(0.95*minvalRat);
961  hUpp->SetTitle(fullTitleRat.c_str());
962 
963  CenterTitles(hDwn);
964  hDwn->SetLineColor(kRed);
965  hDwn->SetLineStyle(2);
966  hDwn->SetMaximum(1.05*maxvalRat);
967  hDwn->SetMinimum(0.95*minvalRat);
968  hDwn->SetTitle(fullTitleRat.c_str());
969 
970  TCanvas* c = new TCanvas(name.c_str(), title.c_str(), 800, 800);
971  TPad* pSpecs = new TPad("pSpecs", "", 0., 0.375, 1., 1.);
972  TPad* pRatio = new TPad("pRatio", "", 0., 0., 1., 0.375);
973  pSpecs->Draw();
974  pRatio->Draw();
975 
976  pSpecs->cd();
977  hErr->Draw("e2");
978  h->Draw("hist same");
979  Simulation();
980 
981  pRatio->cd();
982  hRat->Draw("hist");
983  hUpp->Draw("hist same");
984  hDwn->Draw("hist same");
985 
986  c->Update();
987  out->WriteTObject(c);
988  return;
989  }
990 
991  //----------------------------------------------------------------------------
992  /// Save all of the objects in the input maps to the out directory/file
993  void SaveMaps(
994  TDirectory* out,
995  std::map<std::string, IDecomp*> decomps_nominal,
996  std::map<std::string, std::map<std::string, std::map<int, IDecomp*> > > decomps_shifted,
997  std::map<std::string, PredictionNoExtrap*> predsNE_nominal,
998  std::map<std::string, std::map<std::string, std::map<int, PredictionNoExtrap*> > > predsNE_shifted,
999  std::map<std::string, PredictionSterile*> predsSt_nominal,
1000  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > > predsSt_shifted
1001  )
1002  {
1003  std::map<std::string, TDirectory*> savedirs;
1004  std::string dir;
1005  std::string sep = "__";
1006 
1007  for(const auto& decomp : decomps_nominal) {
1008  dir = "decomp" + sep + decomp.first; // sample name
1009  savedirs[dir] = out->mkdir(dir.c_str());
1010  decomp.second->SaveTo(savedirs[dir]);
1011  }
1012 
1013  for(const auto& sample : decomps_shifted) {
1014  for(const auto& shiftlabel : sample.second) {
1015  for(const auto& decomp : shiftlabel.second) {
1016  dir = "decomp" + sep + sample.first +
1017  sep + shiftlabel.first + sep +
1018  StringFromInt(decomp.first, false);
1019  savedirs[dir] = out->mkdir(dir.c_str());
1020  decomp.second->SaveTo(savedirs[dir]);
1021  }
1022  }
1023  }
1024 
1025  for(const auto& pred : predsNE_nominal) {
1026  dir = "predNE" + sep + pred.first;
1027  savedirs[dir] = out->mkdir(dir.c_str());
1028  pred.second->SaveTo(savedirs[dir]);
1029  }
1030 
1031  for(const auto& sample : predsNE_shifted) {
1032  for(const auto& shiftlabel : sample.second) {
1033  for(const auto& pred : shiftlabel.second) {
1034  dir = "predNE" + sep + sample.first +
1035  sep + shiftlabel.first + sep +
1036  StringFromInt(pred.first, false);
1037  savedirs[dir] = out->mkdir(dir.c_str());
1038  pred.second->SaveTo(savedirs[dir]);
1039  }
1040  }
1041  }
1042 
1043  for(const auto& pred : predsSt_nominal) {
1044  dir = "predSt" + sep + pred.first;
1045  savedirs[dir] = out->mkdir(dir.c_str());
1046  pred.second->SaveTo(savedirs[dir]);
1047  }
1048 
1049  for(const auto& sample : predsSt_shifted) {
1050  for(const auto& shiftlabel : sample.second) {
1051  for(const auto& pred : shiftlabel.second) {
1052  dir = "predSt" + sep + sample.first +
1053  sep + shiftlabel.first + sep +
1054  StringFromInt(pred.first, false);
1055  savedirs[dir] = out->mkdir(dir.c_str());
1056  pred.second->SaveTo(savedirs[dir]);
1057  }
1058  }
1059  }
1060 
1061  return;
1062  }
1063 }
void Simulation()
Definition: tools.h:16
T max(const caf::Proxy< T > &a, T b)
const XML_Char * name
Definition: expat.h:151
int nBins
Definition: plotROC.py:16
std::string fDet
Definition: PPFXHelper.h:104
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::string fComponent
Definition: PPFXHelper.h:103
void SqrtError(TH1 *nom)
Set the errors of the input histogram to be their square roots.
Definition: PPFXHelper.h:139
T sqrt(T number)
Definition: d0nt_math.hpp:156
void SaveMaps(TDirectory *out, std::map< std::string, IDecomp * > decomps_nominal, std::map< std::string, std::map< std::string, std::map< int, IDecomp * > > > decomps_shifted, std::map< std::string, PredictionNoExtrap * > predsNE_nominal, std::map< std::string, std::map< std::string, std::map< int, PredictionNoExtrap * > > > predsNE_shifted, std::map< std::string, PredictionSterile * > predsSt_nominal, std::map< std::string, std::map< std::string, std::map< int, PredictionSterile * > > > predsSt_shifted)
Save all of the objects in the input maps to the out directory/file.
Definition: PPFXHelper.h:1077
static std::unique_ptr< PredictionSterile > LoadFrom(TDirectory *dir, const std::string &name)
void PlotSyst(TH1 *h, TH1 *hp1, TH1 *hm1, TH1 *hp2, TH1 *hm2, TDirectory *out, FILE *text, strings strs)
Definition: PPFXHelper.h:662
void ZeroError(TH1 *nom)
Set the errors of the input histogram to 0.
Definition: PPFXHelper.h:149
Float_t tmp
Definition: plot.C:36
void CenterTitles(TH1 *histo)
Definition: Plots.cxx:1483
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
float abs(float number)
Definition: d0nt_math.hpp:39
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
TH1 * GetNC(IDecomp *specs, double POT)
Definition: PPFXHelper.h:211
const XML_Char * s
Definition: expat.h:262
Charged-current interactions.
Definition: IPrediction.h:39
void AddErrorInQuadrature(TH1 *nom, TH1 *up, bool use_max)
Definition: PPFXHelper.h:163
void prev()
Definition: show_event.C:91
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
void InitializeSystText(FILE *text, strings strs)
Print some initial text about a systematic – the systematic type and sample.
Definition: PPFXHelper.h:263
std::string fSystL
Definition: PPFXHelper.h:109
#define pot
static std::unique_ptr< PredictionExtrap > LoadFrom(TDirectory *dir, const std::string &name)
std::string fPOT
Definition: PPFXHelper.h:105
std::vector< float > Spectrum
Definition: Constants.h:728
double sigma(TH1F *hist, double percentile)
const char sep
std::string StringFromInt(int i, bool withE20)
Helper function that takes an integer and turns it into a string.
Definition: PPFXHelper.h:115
std::vector< double > POT
void PlotSystBand(TH1 *h, TDirectory *out, FILE *text, strings strs)
Definition: PPFXHelper.h:968
std::string fXLabel
Definition: PPFXHelper.h:110
TDirectory * dir
Definition: macro.C:5
std::string fSystType
Definition: PPFXHelper.h:107
Standard interface to all decomposition techniques.
Definition: IDecomp.h:13
Neutral-current interactions.
Definition: IPrediction.h:40
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
A prediction object compatible with sterile oscillations.
std::string fSystS
Definition: PPFXHelper.h:108
void LoadMaps(TDirectory *dir, std::map< std::string, IDecomp * > *nominal, std::map< std::string, std::map< std::string, std::map< int, IDecomp * > > > *shifted)
Definition: PPFXHelper.h:273
const Color_t kTotalMCColor
Definition: Style.h:16
Prediction that just uses FD MC, with no extrapolation.
T min(const caf::Proxy< T > &a, T b)
All neutrinos, any flavor.
Definition: IPrediction.h:26
std::unique_ptr< IDecomp > LoadFrom< IDecomp >(TDirectory *dir, const std::string &label)
Definition: IDecomp.cxx:53
TH1 * GetBG(IDecomp *specs, double POT)
Definition: PPFXHelper.h:236
enum BeamMode kBlue
std::string StringFromDouble(double pot)
Definition: PPFXHelper.h:128
A helper structure to contain a group of string for plotting.
Definition: PPFXHelper.h:101
std::string fSample
Definition: PPFXHelper.h:106
enum BeamMode string