NuSSystFunctions18.h
Go to the documentation of this file.
1 /*
2  This file contains function systematic studies
3  In this file:
4 
5 
6  void SqrtError(TH1* nom);
7  void ZeroError(TH1* nom);
8 
9  void AddErrorInQuadrature(TH1* nom, TH1* up, bool use_max);
10  void AddErrorInQuadrature(TH1* nom, TH1* up, TH1* down, bool use_max);
11  void AddErrorInQuadrature(TH1* nom, TH1* up1, TH1* down1, TH1* up2, TH1* down2, bool use_max);
12 
13  TH1* GetNC(IDecomp* specs, double POT);
14  TH1* GetNC(IPrediction* specs, double POT);
15 
16  TH1* GetBG(IDecomp* specs, double POT);
17  TH1* GetBG(IPrediction* specs, double POT);
18 
19  void InitializeSystText(FILE* text, strings strs);
20 
21  void LoadMaps(TDirectory* dir,
22  std::map<std::string, IDecomp*> * nominal,
23  std::map<std::string, std::map<std::string, std::map<int, IDecomp*> > >* shifted);
24  void LoadMaps(TDirectory* dir,
25  std::map<std::string, PredictionNoExtrap*> * nominal,
26  std::map<std::string, std::map<std::string, std::map<int, PredictionNoExtrap*> > >* shifted);
27  void LoadMaps(TDirectory* dir,
28  std::map<std::string, PredictionSterile*> * nominal,
29  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > >* shifted);
30 
31  void PlotSyst(TH1* h0, TH1* hp1,
32  TDirectory* out, FILE* text, strings strs);
33  void PlotSyst(TH1* h0, TH1* hp1, TH1* hm1, TH1* hp2, TH1* hm2,
34  TDirectory* out, FILE* text, strings strs);
35 
36  void PlotSyst(TH1* nom, IDecomp* shifts,
37  TDirectory* out, FILE* text, strings strs,
38  double POT, bool NC = true, bool use_max);
39  void PlotSyst(TH1* nom, PredictionNoExtrap* shifts,
40  TDirectory* out, FILE* text, strings strs,
41  double POT, bool NC = true, bool use_max);
42  void PlotSyst(TH1* nom, PredictionSterile* shifts,
43  TDirectory* out, FILE* text, strings strs,
44  double POT, bool NC = true, bool use_max);
45 
46  void PlotSyst(TH1* nom, std::map<int, IDecomp*> shifts,
47  TDirectory* out, FILE* text, strings strs,
48  double POT, bool NC = true, bool use_max);
49  void PlotSyst(TH1* nom, std::map<int, PredictionNoExtrap*> shifts,
50  TDirectory* out, FILE* text, strings strs,
51  double POT, bool NC = true, bool use_max);
52  void PlotSyst(TH1* nom, std::map<int, PredictionSterile*> shifts,
53  TDirectory* out, FILE* text, strings strs,
54  double POT, bool NC = true, bool use_max);
55 
56  void PlotSystBand(TH1* h, TDirectory* out, FILE* text, strings strs);
57 
58  void SaveMaps(
59  TDirectory* out,
60  std::map<std::string, IDecomp*> decomps_nominal,
61  std::map<std::string, std::map<std::string, std::map<int, IDecomp*> > > decomps_shifted,
62  std::map<std::string, PredictionNoExtrap*> predsNE_nominal,
63  std::map<std::string, std::map<std::string, std::map<int, PredictionNoExtrap*> > > predsNE_shifted,
64  std::map<std::string, PredictionSterile*> predsSt_nominal,
65  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > > predsSt_shifted
66  );
67 */
68 
69 #pragma once
70 
71 #include <map>
72 #include <iostream>
73 #include <string>
74 #include <vector>
75 
76 #include "CAFAna/Analysis/Calcs.h"
77 #include "CAFAna/Analysis/Style.h"
79 #include "CAFAna/Decomp/IDecomp.h"
85 
88 
91 
92 #include "TCanvas.h"
93 #include "TDirectory.h"
94 #include "TH1.h"
95 #include "TKey.h"
96 #include "TLatex.h"
97 #include "TLegend.h"
98 #include "TString.h"
99 #include "TObjArray.h"
100 #include "TObjString.h"
101 
102 namespace ana
103 {
104  /// A helper structure to contain a group of string for plotting
105  struct strings
106  {
115  };
116 
117  //----------------------------------------------------------------------------
118  /// Set the errors of the input histogram to be their square roots
119  void SqrtError(TH1* nom) {
120  for (int i = 0, n = nom->GetNbinsX(); i <= n; ++i)
121  nom->SetBinError(i, sqrt(nom->GetBinError(i)));
122  return;
123  }
124 
125  //----------------------------------------------------------------------------
126  /// Set the errors of the input histogram to 0
127  void ZeroError(TH1* nom) {
128  for (int i = 0, n = nom->GetNbinsX(); i <= n; ++i)
129  nom->SetBinError(i, 0.);
130  return;
131  }
132 
133  //----------------------------------------------------------------------------
134  /// Combine errors in quadrature, setting error bins in nominal spectra
135  /// Leaves errors squared -- square root must be applied later
136  /// \param nom The nominal histogram
137  /// \param up The shifted spectra
138  /// \param use_max Only keep the max error, don't combine in quadrature
139  void AddErrorInQuadrature(TH1* nom, TH1* up, bool use_max) {
140 
141  double error;
142  for (int i = 0, n = nom->GetNbinsX(); i <= n; ++i) {
143  error = std::abs(up->GetBinContent(i) - nom->GetBinContent(i));
144  double error_sq = error*error;
145  double err_prev = std::abs(nom->GetBinError(i));
146  if (use_max)
147  nom->SetBinError(i, std::max(err_prev, error_sq));
148  else
149  nom->SetBinError(i, err_prev + error_sq);
150  }
151 
152  return;
153  }
154 
155  //----------------------------------------------------------------------------
156  /// Combine errors in quadrature, setting error bins in nominal spectra
157  /// Takes the larger difference from nominal from either up or down shift
158  /// Leaves errors squared -- square root must be applied later
159  /// \param nom The nominal histogram
160  /// \param up The +1 sigma shifted spectra
161  /// \param down The -1 sigma shifted spectra
162  /// \param use_max Only keep the max error, don't combine in quadrature
163  void AddErrorInQuadrature(TH1* nom, TH1* up, TH1* down, bool use_max) {
164 
165  double error;
166  for (int i = 0, n = nom->GetNbinsX(); i <= n; ++i) {
167  error = std::max(std::abs(up ->GetBinContent(i) - nom ->GetBinContent(i)),
168  std::abs(nom->GetBinContent(i) - down->GetBinContent(i)));
169  double error_sq = error*error;
170  double err_prev = std::abs(nom->GetBinError(i));
171  if (use_max)
172  nom->SetBinError(i, std::max(err_prev, error_sq));
173  else
174  nom->SetBinError(i, err_prev + error_sq);
175  }
176 
177  return;
178  }
179 
180  //----------------------------------------------------------------------------
181  /// Combine errors in quadrature, setting error bins in nominal spectra
182  /// Takes the larger difference from nominal from either up or down shift
183  /// Leaves errors squared -- square root must be applied later
184  /// This version takes multiple up/down hists (for example, when splitting
185  /// up systematics)
186  /// \param nom The nominal histogram
187  /// \param up1/up2 The +1 sigma shifted spectra
188  /// \param down1/down2 The -1 sigma shifted spectra
189  /// \param use_max Only keep the max error, don't combine in quadrature
190  void AddErrorInQuadrature(TH1* nom, TH1* up1, TH1* down1, TH1* up2, TH1* down2, bool use_max = false) {
191 
192  double error1, error2;
193  for(int i = 0, n = nom->GetNbinsX(); i <= n; ++i) {
194 
195  error1 = std::max(std::abs(up1 ->GetBinContent(i) - nom ->GetBinContent(i)),
196  std::abs(nom->GetBinContent(i) - down1->GetBinContent(i)));
197  double error1_sq = error1*error1;
198 
199  error2 = std::max(std::abs(up2 ->GetBinContent(i) - nom ->GetBinContent(i)),
200  std::abs(nom->GetBinContent(i) - down2->GetBinContent(i)));
201  double error2_sq = error2*error2;
202 
203  double err_prev = std::abs(nom->GetBinError(i));
204  if(use_max)
205  nom->SetBinError(i, std::max(err_prev, error1_sq));
206  else
207  nom->SetBinError(i, err_prev + error1_sq + error2_sq);
208 
209  }
210 
211  return;
212  }
213 
214  //----------------------------------------------------------------------------
215  /// Get the NC spectrum from a prediction object at a specific POT value
216  /// Sets errors to 0
217  TH1* GetNC(IPrediction* specs, double POT) {
218 
219  //osc::OscCalculatorPMNSOpt* calc = (osc::OscCalculatorPMNSOpt*)DefaultOscCalc();
220  //osc::OscCalculatorSterile* calc = DefaultSterileCalc(3);
222 
224 
225  TH1* ret = s.ToTH1(POT);
226  ZeroError(ret);
227 
228  return ret;
229 
230  }
231 
232  //----------------------------------------------------------------------------
233  /// Get the background spectrum from a prediction object at a specific POT value
234  /// Sets errors to 0
235  TH1* GetBG(IPrediction* specs, double POT) {
236 
237  //osc::OscCalculatorPMNSOpt* calc = (osc::OscCalculatorPMNSOpt*)DefaultOscCalc();
239 
241 
242  TH1* ret = s.ToTH1(POT);
243  ZeroError(ret);
244 
245  return ret;
246 
247  }
248 
249  //----------------------------------------------------------------------------
250  /// Print some initial text about a systematic -- the systematic type and sample
251  void InitializeSystText(FILE* text, strings strs) {
252  fprintf(text, "%s %s\n\n", strs.fSystType.c_str(), strs.fSample.c_str());
253  return;
254  }
255 
256  //----------------------------------------------------------------------------
257  /// Load nominal and shifted NDPredictionSterile objects from a TDirectory
258  /// The maps to contain these objects must already exist
259  void LoadMaps(TDirectory* dir,
260  std::map<std::string, NDPredictionSterile*>* nominal,
261  std::map<std::string, std::map<std::string, std::map<int, NDPredictionSterile*> > >* shifted) {
262 
263  TDirectory* tmp = gDirectory;
264  TDirectory* keyDir = gDirectory;
265  dir->cd();
266  TIter nextkey(dir->GetListOfKeys());
267  TKey* key;
268  TKey* oldkey = 0;
269 
270  std::string sep = "__";
271  std::string cafanaType = "NDPred";
272  std::string sampleLabel = "";
273  std::string systLabel = "";
274  int sigma = 0;
275 
276  // Loop through the contents of the TDirectory
277  while ((key = (TKey*)nextkey())) {
278  if (oldkey && !strcmp(oldkey->GetName(), key->GetName()))
279  continue;
280  if (key->ReadObj()->IsFolder())
281  keyDir = (TDirectory*)key->ReadObj();
282  else
283  continue;
284 
285  // tokenize!
286  std::string name(key->GetName());
287  std::vector<std::string> tokens;
288  std::size_t pos = 0, prev = 0;
289  while (pos != std::string::npos) { // Loop over the TObject name
290  pos = name.find(sep, prev); // Find the next separator
291  if (pos != std::string::npos) {
292  tokens.push_back(name.substr(prev, pos-prev));
293  prev = pos+2; // Move past the current separator
294  }
295  else
296  tokens.push_back(name.substr(prev, name.length()-prev));
297  }
298 
299  // Make sure the object we found is the type we want
300  if (cafanaType.compare(tokens[0]) != 0)
301  continue;
302 
303  // Fill the correct elements of the maps with the prediction object
304  sampleLabel = tokens[1];
305  if (tokens.size() > 2) {
306  systLabel = tokens[2];
307  sigma = std::stoi(tokens[3]);
308  (*shifted)[sampleLabel][systLabel][sigma] = NDPredictionSterile::LoadFrom(keyDir).release();
309  }
310  else
311  (*nominal)[sampleLabel] = NDPredictionSterile::LoadFrom(keyDir).release();
312 
313  }
314 
315  tmp->cd();
316  return;
317 
318  }
319 
320  //----------------------------------------------------------------------------
321  /// Load nominal and shifted FDPredictionSterile objects from a TDirectory
322  /// The maps to contain these objects must already exist
323  void LoadMaps(TDirectory* dir,
324  std::map<std::string, FDPredictionSterile*>* nominal,
325  std::map<std::string, std::map<std::string, std::map<int, FDPredictionSterile*> > >* shifted) {
326 
327  TDirectory* tmp = gDirectory;
328  TDirectory* keyDir = gDirectory;
329  dir->cd();
330  TIter nextkey(dir->GetListOfKeys());
331  TKey* key;
332  TKey* oldkey = 0;
333 
334  std::string sep = "__";
335  std::string cafanaType = "FDPred";
336  std::string sampleLabel = "";
337  std::string systLabel = "";
338  int sigma = 0;
339 
340  // Loop through the contents of the TDirectory
341  while ((key = (TKey*)nextkey())) {
342  if (oldkey && !strcmp(oldkey->GetName(), key->GetName()))
343  continue;
344  if (key->ReadObj()->IsFolder())
345  keyDir = (TDirectory*)key->ReadObj();
346  else
347  continue;
348 
349  // tokenize!
350  std::string name(key->GetName());
351  std::vector<std::string> tokens;
352  std::size_t pos = 0, prev = 0;
353  while (pos != std::string::npos) { // Loop over the TObject name
354  pos = name.find(sep, prev); // Find the next separator
355  if (pos != std::string::npos) {
356  tokens.push_back(name.substr(prev, pos-prev));
357  prev = pos+2; // Move past the current separator
358  }
359  else
360  tokens.push_back(name.substr(prev, name.length()-prev));
361  }
362 
363  // Make sure the object we found is the type we want
364  if (cafanaType.compare(tokens[0]) != 0)
365  continue;
366 
367  // Fill the correct elements of the maps with the prediction object
368  sampleLabel = tokens[1];
369  if (tokens.size() > 2) {
370  systLabel = tokens[2];
371  sigma = std::stoi(tokens[3]);
372  (*shifted)[sampleLabel][systLabel][sigma] = FDPredictionSterile::LoadFrom(keyDir).release();
373  }
374  else
375  (*nominal)[sampleLabel] = FDPredictionSterile::LoadFrom(keyDir).release();
376 
377  }
378 
379  tmp->cd();
380  return;
381 
382  }
383 
384  //----------------------------------------------------------------------------
385  /// Load nominal and shifted PredictionSterile objects from a TDirectory
386  /// The maps to contain these objects must already exist
387  void LoadMaps(TDirectory* dir,
388  std::map<std::string, PredictionSterile*>* nominal,
389  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > >* shifted) {
390 
391  TDirectory* tmp = gDirectory;
392  TDirectory* keyDir = gDirectory;
393  dir->cd();
394  TIter nextkey(dir->GetListOfKeys());
395  TKey* key;
396  TKey* oldkey = 0;
397 
398  std::string sep = "__";
399  std::string cafanaType = "ExtrapPred";
400  std::string sampleLabel = "";
401  std::string systLabel = "";
402  int sigma = 0;
403 
404  // Loop through the contents of the TDirectory
405  while ((key = (TKey*)nextkey())) {
406  if (oldkey && !strcmp(oldkey->GetName(), key->GetName()))
407  continue;
408  if (key->ReadObj()->IsFolder())
409  keyDir = (TDirectory*)key->ReadObj();
410  else
411  continue;
412 
413  // tokenize!
414  std::string name(key->GetName());
415  std::vector<std::string> tokens;
416  std::size_t pos = 0, prev = 0;
417  while (pos != std::string::npos) { // Loop over the TObject name
418  pos = name.find(sep, prev); // Find the next separator
419  if (pos != std::string::npos) {
420  tokens.push_back(name.substr(prev, pos-prev));
421  prev = pos+2; // Move past the current separator
422  }
423  else
424  tokens.push_back(name.substr(prev, name.length()-prev));
425  }
426 
427  // Make sure the object we found is the type we want
428  if (cafanaType.compare(tokens[0]) != 0)
429  continue;
430 
431  // Fill the correct elements of the maps with the prediction object
432  sampleLabel = tokens[1];
433  if (tokens.size() > 2) {
434  systLabel = tokens[2];
435  sigma = std::stoi(tokens[3]);
436  (*shifted)[sampleLabel][systLabel][sigma] = PredictionSterile::LoadFrom(keyDir).release();
437  }
438  else
439  (*nominal)[sampleLabel] = PredictionSterile::LoadFrom(keyDir).release();
440 
441  }
442 
443  tmp->cd();
444  return;
445 
446  }
447 
448  //----------------------------------------------------------------------------
449  /// Plot a nominal spectrum and a shifted spectra
450  /// This is for an on/off sytle systematic
451  /// \param h The nominal spectrum (systematic off)
452  /// \param hp1 The shifted spectrum (systematic on)
453  /// \param out The directory/file for root output
454  /// \param text The file for text output
455  /// \param strs Contains strings for pretty plots
456  void PlotSyst(TH1* h, TH1* hp1, TDirectory* out, FILE* text, strings strs)
457  {
458  // Construct histogram names and titles
459  std::string xlabel = strs.fXLabel;
460  std::string ylabel = "Events / " + strs.fPOT + " / 0.25 GeV";
461  std::string title = strs.fDet + " " + strs.fComponent +
462  " Error for " + strs.fSystL + " Systematic";
463  std::string fullTitle = ";" + xlabel + ";" + ylabel;
464  std::string name = "c" + strs.fComponent + strs.fDet +
465  strs.fSystS + strs.fSample;
466  std::string ylabelRat = "Shifted / Nominal";
467  std::string fullTitleRat = ";" + xlabel + ";" + ylabelRat;
468 
469  TH1* hRat = (TH1*)h->Clone();
470  TH1* hp1Rat = (TH1*)hp1->Clone();
471 
472  // Determine overall percentage difference
473  // Rate only, no shape
474  int nBins = h->GetNbinsX();
475  double percent_diff = 100. * std::abs(hp1->Integral(1, nBins) - h->Integral(1, nBins))/
476  h->Integral(1, nBins);
477 
478  // Write to file
479  fprintf(text, "%s %s %.2s: +/-%6.4f%%\n",
480  strs.fSystS.c_str(), strs.fComponent.c_str(), strs.fDet.c_str(),
481  percent_diff);
482 
483  // Get ratios of nominal (1. in each bin) and
484  // shifted (divded by nominal)
485  for(int i = 0; i <= nBins+1; ++i)
486  hRat->SetBinContent(i, 1.);
487  hp1Rat->Divide(h);
488 
489  // Get max and min values
490  double maxval = h->GetBinContent(h->GetMaximumBin());
491  double minval = h->GetBinContent(h->GetMinimumBin());
492  maxval = std::max(maxval, hp1->GetBinContent(hp1->GetMaximumBin()));
493  minval = std::min(minval, hp1->GetBinContent(hp1->GetMinimumBin()));
494  if (maxval < minval)
495  std::swap(maxval, minval);
496 
497  // Get max and min ratio values
498  double maxvalRat = hRat->GetBinContent(hRat->GetMaximumBin());
499  double minvalRat = hRat->GetBinContent(hRat->GetMinimumBin());
500  maxvalRat = std::max(maxvalRat, hp1Rat->GetBinContent(hp1Rat->GetMaximumBin()));
501  for(int i = 1; i <= nBins; ++i)
502  if (hp1Rat->GetBinContent(i) > 0.)
503  minvalRat = std::min(minvalRat, hp1Rat->GetBinContent(i));
504  if (std::abs(maxvalRat - 1.) > std::abs(1. - minvalRat)) {
505  minvalRat = 1 - std::abs(maxvalRat - 1.);
506  maxvalRat = 1 + std::abs(maxvalRat - 1.);
507  } else {
508  minvalRat = 1 - std::abs(1. - minvalRat);
509  maxvalRat = 1 + std::abs(1. - minvalRat);
510  }
511 
512  // Formatting
513  CenterTitles(h);
514  h->SetLineColor(kBlack);
515  h->SetMaximum(1.1*maxval);
516  h->SetMinimum(0);
517  h->SetTitle(fullTitle.c_str());
518 
519  CenterTitles(hp1);
520  hp1->SetLineColor(kRed);
521  hp1->SetLineStyle(2);
522  hp1->SetMaximum(1.1*maxval);
523  hp1->SetMinimum(0);
524  hp1->SetTitle(fullTitle.c_str());
525 
526  CenterTitles(hRat);
527  hRat->SetLineColor(kBlack);
528  hRat->SetMaximum(1.1*maxvalRat);
529  hRat->SetMinimum(0.9*minvalRat);
530  hRat->SetTitle(fullTitleRat.c_str());
531 
532  CenterTitles(hp1Rat);
533  hp1Rat->SetLineColor(kRed);
534  hp1Rat->SetLineStyle(2);
535  hp1Rat->SetMaximum(1.1*maxvalRat);
536  hp1Rat->SetMinimum(0.9*minvalRat);
537  hp1Rat->SetTitle(fullTitleRat.c_str());
538 
539  // Legend
540  double xL = 0.6, xR = 0.85;
541  double yB = 0.6, yT = 0.85;
542  TLegend* leg = new TLegend(xL, yB, xR, yT);
543  leg->AddEntry(h, "Nominal", "l");
544  leg->AddEntry(hp1, "Shifted", "l");
545 
546  // Draw all hists
547  TCanvas* c = new TCanvas(name.c_str(), title.c_str(), 800, 800);
548  TPad* pSpecs = new TPad("pSpecs", "", 0., 0.375, 1., 1.);
549  TPad* pRatio = new TPad("pRatio", "", 0., 0., 1., 0.375);
550  pSpecs->Draw();
551  pRatio->Draw();
552 
553  pSpecs->cd();
554  h->Draw("hist");
555  hp1->Draw("hist same");
556  leg->Draw();
557  Simulation();
558 
559  pRatio->cd();
560  hRat->Draw("hist");
561  hp1Rat->Draw("hist same");
562 
563  c->Update();
564  out->WriteTObject(c);
565  return;
566 
567  }
568 
569  //----------------------------------------------------------------------------
570  /// Plot a nominal spectrum and a shifted spectra
571  /// This is for a systematic with multiple sigma levels
572  /// \param h The nominal spectrum
573  /// \param hp1 The shifted spectrum, +1 sigma
574  /// \param hm1 The shifted spectrum, -1 sigma
575  /// \param hp2 The shifted spectrum, +2 sigma
576  /// \param hm2 The shifted spectrum, -2 sigma
577  /// \param out The directory/file for root output
578  /// \param text The file for text output
579  /// \param strs Contains strings for pretty plots
580  void PlotSyst(TH1* h, TH1* hp1, TH1* hm1, TH1* hp2, TH1* hm2,
581  TDirectory* out, FILE* text, strings strs)
582  {
583  std::string xlabel = strs.fXLabel;
584  std::string ylabel = "Events / " + strs.fPOT + " / 0.25 GeV";
585  std::string title = strs.fDet + " " + strs.fComponent +
586  " Error for " + strs.fSystL + " Systematic";
587  std::string fullTitle = ";" + xlabel + ";" + ylabel;
588  std::string name = "c" + strs.fComponent + strs.fDet +
589  strs.fSystS + strs.fSample;
590  std::string ylabelRat = "Shifted / Nominal";
591  std::string fullTitleRat = ";" + xlabel + ";" + ylabelRat;
592 
593  // Ratio hists
594  TH1* hRat = (TH1*)h->Clone();
595  TH1* hp1Rat = (TH1*)hp1->Clone();
596  TH1* hm1Rat = (TH1*)hm1->Clone();
597  TH1* hp2Rat = (TH1*)hp2->Clone();
598  TH1* hm2Rat = (TH1*)hm2->Clone();
599 
600  // Get absolute percentage differences between the spectra
601  // No shape info
602  // Just +/- 1sigma
603  int nBins = h->GetNbinsX();
604  double nom_evts = h->Integral(1, nBins);
605  double evts_diff = 0.;
606  double percent_diff = 100.;
607  double percent_diff_p = 100.;
608  double percent_diff_m = 100.;
609  for (int bin = 0; bin <= nBins+1; ++bin) {
610  hRat->SetBinContent(bin, 1.);
611  if (bin != 0 && bin != nBins+1)
612  evts_diff += std::max(std::abs(hp1->GetBinContent(bin) - h->GetBinContent(bin)),
613  std::abs(hm1->GetBinContent(bin) - h->GetBinContent(bin)));
614  }
615  percent_diff *= (evts_diff/nom_evts);
616  percent_diff_p *= ((hp1->Integral(1, nBins) - nom_evts)/nom_evts);
617  percent_diff_m *= ((nom_evts - hm1->Integral(1, nBins))/nom_evts);
618 
619  fprintf(text, "%s %.2s %.2s: +/-%6.4f%%, +%6.4f%%, -%6.4f%%\n",
620  strs.fSystS.c_str(), strs.fComponent.c_str(), strs.fDet.c_str(),
621  percent_diff, percent_diff_p, percent_diff_m);
622 
623  // Get ratios
624  hp1Rat->Divide(h);
625  hm1Rat->Divide(h);
626  hp2Rat->Divide(h);
627  hm2Rat->Divide(h);
628 
629  // Get max and min values
630  double maxval = h->GetBinContent(h->GetMaximumBin());
631  maxval = std::max(maxval, hp1->GetBinContent(hp1->GetMaximumBin()));
632  maxval = std::max(maxval, hm1->GetBinContent(hm1->GetMaximumBin()));
633  maxval = std::max(maxval, hp2->GetBinContent(hp2->GetMaximumBin()));
634  maxval = std::max(maxval, hm2->GetBinContent(hm2->GetMaximumBin()));
635  double minval = h->GetBinContent(h->GetMinimumBin());
636  minval = std::min(minval, hp1->GetBinContent(hp1->GetMinimumBin()));
637  minval = std::min(minval, hm1->GetBinContent(hm1->GetMinimumBin()));
638  minval = std::min(minval, hp2->GetBinContent(hp2->GetMinimumBin()));
639  minval = std::min(minval, hm2->GetBinContent(hm2->GetMinimumBin()));
640  if (maxval < minval)
641  std::swap(maxval, minval);
642 
643  double maxvalRat = hRat->GetBinContent(hRat->GetMaximumBin());
644  maxvalRat = std::max(maxvalRat, hp1Rat->GetBinContent(hp1Rat->GetMaximumBin()));
645  maxvalRat = std::max(maxvalRat, hm1Rat->GetBinContent(hm1Rat->GetMaximumBin()));
646  maxvalRat = std::max(maxvalRat, hp2Rat->GetBinContent(hp2Rat->GetMaximumBin()));
647  maxvalRat = std::max(maxvalRat, hm2Rat->GetBinContent(hm2Rat->GetMaximumBin()));
648  double minvalRat = hRat->GetBinContent(hRat->GetMinimumBin());
649  for (int i = 1; i <= nBins; ++i) {
650  if (hp1Rat->GetBinContent(i) > 0.)
651  minvalRat = std::min(minvalRat, hp1Rat->GetBinContent(i));
652  if (hm1Rat->GetBinContent(i) > 0.)
653  minvalRat = std::min(minvalRat, hm1Rat->GetBinContent(i));
654  if (hp2Rat->GetBinContent(i) > 0.)
655  minvalRat = std::min(minvalRat, hp2Rat->GetBinContent(i));
656  if (hm2Rat->GetBinContent(i) > 0.)
657  minvalRat = std::min(minvalRat, hm2Rat->GetBinContent(i));
658  }
659  if(std::abs(maxvalRat - 1.) > std::abs(1. - minvalRat)) {
660  minvalRat = 1 - std::abs(maxvalRat - 1.);
661  maxvalRat = 1 + std::abs(maxvalRat - 1.);
662  } else {
663  minvalRat = 1 - std::abs(1. - minvalRat);
664  maxvalRat = 1 + std::abs(1. - minvalRat);
665  }
666 
667  // Draw
668  CenterTitles(h);
669  h->SetLineColor(kBlack);
670  h->SetMaximum(1.1*maxval);
671  h->SetMinimum(0);
672  h->SetTitle(fullTitle.c_str());
673 
674  CenterTitles(hp1);
675  hp1->SetLineColor(kRed);
676  hp1->SetLineStyle(2);
677  hp1->SetMaximum(1.1*maxval);
678  hp1->SetMinimum(0);
679  hp1->SetTitle(fullTitle.c_str());
680 
681  CenterTitles(hm1);
682  hm1->SetLineColor(kRed);
683  hm1->SetLineStyle(2);
684  hm1->SetMaximum(1.1*maxval);
685  hm1->SetMinimum(0);
686  hm1->SetTitle(fullTitle.c_str());
687 
688  CenterTitles(hp2);
689  hp2->SetLineColor(kBlue);
690  hp2->SetLineStyle(2);
691  hp2->SetMaximum(1.1*maxval);
692  hp2->SetMinimum(0);
693  hp2->SetTitle(fullTitle.c_str());
694 
695  CenterTitles(hm2);
696  hm2->SetLineColor(kBlue);
697  hm2->SetLineStyle(2);
698  hm2->SetMaximum(1.1*maxval);
699  hm2->SetMinimum(0);
700  hm2->SetTitle(fullTitle.c_str());
701 
702  CenterTitles(hRat);
703  hRat->SetLineColor(kBlack);
704  hRat->SetMaximum(1.1*maxvalRat);
705  hRat->SetMinimum(0.9*minvalRat);
706  hRat->SetTitle(fullTitleRat.c_str());
707 
708  CenterTitles(hp1Rat);
709  hp1Rat->SetLineColor(kRed);
710  hp1Rat->SetLineStyle(2);
711  hp1Rat->SetMaximum(1.1*maxvalRat);
712  hp1Rat->SetMinimum(0.9*minvalRat);
713  hp1Rat->SetTitle(fullTitleRat.c_str());
714 
715  CenterTitles(hm1Rat);
716  hm1Rat->SetLineColor(kRed);
717  hm1Rat->SetLineStyle(2);
718  hm1Rat->SetMaximum(1.1*maxvalRat);
719  hm1Rat->SetMinimum(0.9*minvalRat);
720  hm1Rat->SetTitle(fullTitleRat.c_str());
721 
722  CenterTitles(hp2Rat);
723  hp2Rat->SetLineColor(kBlue);
724  hp2Rat->SetLineStyle(2);
725  hp2Rat->SetMaximum(1.1*maxvalRat);
726  hp2Rat->SetMinimum(0.9*minvalRat);
727  hp2Rat->SetTitle(fullTitleRat.c_str());
728 
729  CenterTitles(hm2Rat);
730  hm2Rat->SetLineColor(kBlue);
731  hm2Rat->SetLineStyle(2);
732  hm2Rat->SetMaximum(1.1*maxvalRat);
733  hm2Rat->SetMinimum(0.9*minvalRat);
734  hm2Rat->SetTitle(fullTitleRat.c_str());
735 
736  double xL = 0.6, xR = 0.85;
737  double yB = 0.6, yT = 0.85;
738  TLegend* leg = new TLegend(xL, yB, xR, yT);
739  leg->AddEntry(h, "Nominal", "l");
740  leg->AddEntry(hp1, "+/- 1 #sigma", "l");
741  leg->AddEntry(hp2, "+/- 2 #sigma", "l");
742 
743  TCanvas* c = new TCanvas(name.c_str(), title.c_str(), 800, 800);
744  TPad* pSpecs = new TPad("pSpecs", "", 0., 0.375, 1., 1.);
745  TPad* pRatio = new TPad("pRatio", "", 0., 0., 1., 0.375);
746  pSpecs->Draw();
747  pRatio->Draw();
748 
749  pSpecs->cd();
750  h->Draw("hist");
751  hp2->Draw("hist same");
752  hm2->Draw("hist same");
753  hp1->Draw("hist same");
754  hm1->Draw("hist same");
755  leg->Draw();
756  Simulation();
757 
758  pRatio->cd();
759  hRat->Draw("hist");
760  hp2Rat->Draw("hist same");
761  hm2Rat->Draw("hist same");
762  hp1Rat->Draw("hist same");
763  hm1Rat->Draw("hist same");
764  c->Update();
765  out->WriteTObject(c);
766 
767  return;
768 
769  }
770 
771  //----------------------------------------------------------------------------
772  /// Plot a nominal spectrum and a shifted spectra
773  /// This is for a systematic with multiple sigma levels
774  /// \param h The nominal spectrum
775  /// \param hp1 The shifted spectrum, +1 sigma
776  /// \param hm1 The shifted spectrum, -1 sigma
777  /// \param out The directory/file for root output
778  /// \param text The file for text output
779  /// \param strs Contains strings for pretty plots
780  void PlotSyst(TH1* h, TH1* hp1, TH1* hm1,
781  TDirectory* out, FILE* text, strings strs) {
782 
783  std::string xlabel = strs.fXLabel;
784  std::string ylabel = "Events / " + strs.fPOT + " / 0.25 GeV";
785  std::string title = strs.fDet + " " + strs.fComponent +
786  " Error for " + strs.fSystL + " Systematic";
787  std::string fullTitle = ";" + xlabel + ";" + ylabel;
788  std::string name = "c" + strs.fComponent + strs.fDet +
789  strs.fSystS + strs.fSample;
790  std::string ylabelRat = "Shifted / Nominal";
791  std::string fullTitleRat = ";" + xlabel + ";" + ylabelRat;
792 
793  TH1* hRat = (TH1*)h->Clone();
794  TH1* hp1Rat = (TH1*)hp1->Clone();
795  TH1* hm1Rat = (TH1*)hm1->Clone();
796 
797  int nBins = h->GetNbinsX();
798  double nom_evts = h->Integral(1, nBins);
799  double evts_diff = 0.;
800  double percent_diff = 100.;
801  double percent_diff_p = 100.;
802  double percent_diff_m = 100.;
803  for(int i = 0; i <= nBins+1; ++i) {
804  hRat->SetBinContent(i, 1.);
805 
806  if(i != 0 && i != nBins+1) {
807  evts_diff += std::max(std::abs(hp1->GetBinContent(i) - h->GetBinContent(i)),
808  std::abs(hm1->GetBinContent(i) - h->GetBinContent(i)));
809  }
810  }
811  percent_diff *= (evts_diff/nom_evts);
812  percent_diff_p *= ((hp1->Integral(1, nBins) - nom_evts)/nom_evts);
813  percent_diff_m *= ((nom_evts - hm1->Integral(1, nBins))/nom_evts);
814 
815  fprintf(text, "%s %.2s %.2s: +/-%6.4f%%, +%6.4f%%, -%6.4f%%\n",
816  strs.fSystS.c_str(), strs.fComponent.c_str(), strs.fDet.c_str(),
817  percent_diff, percent_diff_p, percent_diff_m);
818 
819  hp1Rat->Divide(h);
820  hm1Rat->Divide(h);
821 
822  double maxval = h->GetBinContent(h->GetMaximumBin());
823  maxval = std::max(maxval, hp1->GetBinContent(hp1->GetMaximumBin()));
824  maxval = std::max(maxval, hm1->GetBinContent(hm1->GetMaximumBin()));
825 
826  double minval = h->GetBinContent(h->GetMinimumBin());
827  minval = std::min(minval, hp1->GetBinContent(hp1->GetMinimumBin()));
828  minval = std::min(minval, hm1->GetBinContent(hm1->GetMinimumBin()));
829 
830  if(maxval < minval) { std::swap(maxval, minval); }
831 
832  double maxvalRat = hRat->GetBinContent(hRat->GetMaximumBin());
833  maxvalRat = std::max(maxvalRat, hp1Rat->GetBinContent(hp1Rat->GetMaximumBin()));
834  maxvalRat = std::max(maxvalRat, hm1Rat->GetBinContent(hm1Rat->GetMaximumBin()));
835 
836  double minvalRat = hRat->GetBinContent(hRat->GetMinimumBin());
837  for(int i = 1; i <= nBins; ++i) {
838  if(hp1Rat->GetBinContent(i) > 0.) {
839  minvalRat = std::min(minvalRat, hp1Rat->GetBinContent(i));
840  }
841  if(hm1Rat->GetBinContent(i) > 0.) {
842  minvalRat = std::min(minvalRat, hm1Rat->GetBinContent(i));
843  }
844  }
845 
846  if(std::abs(maxvalRat - 1.) > std::abs(1. - minvalRat)) {
847  minvalRat = 1 - std::abs(maxvalRat - 1.);
848  maxvalRat = 1 + std::abs(maxvalRat - 1.);
849  }
850  else {
851  minvalRat = 1 - std::abs(1. - minvalRat);
852  maxvalRat = 1 + std::abs(1. - minvalRat);
853  }
854 
855  CenterTitles(h);
856  h->SetLineColor(kBlack);
857  h->SetMaximum(1.1*maxval);
858  h->SetMinimum(0);
859  h->SetTitle(fullTitle.c_str());
860 
861  CenterTitles(hp1);
862  hp1->SetLineColor(kRed);
863  hp1->SetLineStyle(2);
864  hp1->SetMaximum(1.1*maxval);
865  hp1->SetMinimum(0);
866  hp1->SetTitle(fullTitle.c_str());
867 
868  CenterTitles(hm1);
869  hm1->SetLineColor(kRed);
870  hm1->SetLineStyle(2);
871  hm1->SetMaximum(1.1*maxval);
872  hm1->SetMinimum(0);
873  hm1->SetTitle(fullTitle.c_str());
874 
875  CenterTitles(hRat);
876  hRat->SetLineColor(kBlack);
877  hRat->SetMaximum(1.1*maxvalRat);
878  hRat->SetMinimum(0.9*minvalRat);
879  hRat->SetTitle(fullTitleRat.c_str());
880 
881  CenterTitles(hp1Rat);
882  hp1Rat->SetLineColor(kRed);
883  hp1Rat->SetLineStyle(2);
884  hp1Rat->SetMaximum(1.1*maxvalRat);
885  hp1Rat->SetMinimum(0.9*minvalRat);
886  hp1Rat->SetTitle(fullTitleRat.c_str());
887 
888  CenterTitles(hm1Rat);
889  hm1Rat->SetLineColor(kRed);
890  hm1Rat->SetLineStyle(2);
891  hm1Rat->SetMaximum(1.1*maxvalRat);
892  hm1Rat->SetMinimum(0.9*minvalRat);
893  hm1Rat->SetTitle(fullTitleRat.c_str());
894 
895  double xL = 0.6, xR = 0.85;
896  double yB = 0.6, yT = 0.85;
897  TLegend* leg = new TLegend(xL, yB, xR, yT);
898  leg->AddEntry(h, "Nominal", "l");
899  leg->AddEntry(hp1, "+/- 1 #sigma", "l");
900 
901  TCanvas* c = new TCanvas(name.c_str(), title.c_str(), 800, 800);
902  TPad* pSpecs = new TPad("pSpecs", "", 0., 0.375, 1., 1.);
903  TPad* pRatio = new TPad("pRatio", "", 0., 0., 1., 0.375);
904  pSpecs->Draw();
905  pRatio->Draw();
906 
907  pSpecs->cd();
908  h->Draw("hist");
909  hp1->Draw("hist same");
910  hm1->Draw("hist same");
911  leg->Draw();
912  Simulation();
913 
914  pRatio->cd();
915  hRat->Draw("hist");
916  hp1Rat->Draw("hist same");
917  hm1Rat->Draw("hist same");
918  c->Update();
919  out->WriteTObject(c);
920 
921  return;
922 
923  }
924 
925  //----------------------------------------------------------------------------
926  /// Plot a nominal spectrum and a shifted spectra
927  /// This is for a systematic with light level method
928  /// \param h The nominal spectrum
929  /// \param hup The shifted spectrum, +1 sigma
930  /// \param hdn The shifted spectrum, -1 sigma
931  /// \param hsx The shape varied spectrum
932  /// \param out The directory/file for root output
933  /// \param text The file for text output
934  /// \param strs Contains strings for pretty plots
935  void PlotSyst(TH1* nom, TH1* hup, TH1* hdn, TH1* hsx,
936  TDirectory* rootOUT, FILE* textOFS,
937  strings strs) {
938 
939  double nNom = nom->Integral();
940  double n5Up = hup->Integral();
941  double n5Dn = hdn->Integral();
942  double nSlX = hsx->Integral();
943 
944  fprintf(textOFS, "\nEvent counts for %.2s %.2s %s Systs:\n",
945  strs.fDet.c_str(), strs.fComponent.c_str(), strs.fSystType.c_str());
946  fprintf(textOFS, "NC: %6.2f, Flat Up: %6.2f, Flat Down: %6.2f, Slope X: %6.2f",
947  nNom, n5Up, n5Dn, nSlX);
948  fprintf(textOFS, "Flat Up Diff: %6.2f, Flat Down Diff: %6.2f, Slope X Diff: %6.2f",
949  100.*std::abs(n5Up/nNom - 1.), 100.*std::abs(n5Dn/nNom - 1.),
950  100.*std::abs(nSlX/nNom - 1.));
951 
952  std::string xlabel = strs.fXLabel;
953  std::string ylabel = "Events / " + strs.fPOT + " / 0.25 GeV";
954  std::string title = strs.fDet + " " + strs.fComponent +
955  " Errors for All " + strs.fSystType + " Systematics";
956  std::string fullTitle = ";" + xlabel + ";" + ylabel;
957  std::string name = "c" + strs.fComponent + strs.fDet +
958  strs.fSystType + "SystsAll" + strs.fSample;
959  std::string ylabelRat = "Shifted / Nominal";
960  std::string fullTitleRat = ";" + xlabel + ";" + ylabelRat;
961 
962  TH1* hRat = (TH1*)nom->Clone();
963  TH1* hUpRat = (TH1*)hup->Clone();
964  TH1* hDnRat = (TH1*)hdn->Clone();
965  TH1* hSXRat = (TH1*)hsx->Clone();
966 
967  int nBins = nom->GetNbinsX();
968  for (int i = 0; i <= nBins+1; ++i)
969  hRat->SetBinContent(i, 1.);
970 
971  hUpRat->Divide(nom);
972  hDnRat->Divide(nom);
973  hSXRat->Divide(nom);
974 
975  double maxval = nom->GetBinContent(nom->GetMaximumBin());
976  maxval = std::max(maxval, hup->GetBinContent(hup->GetMaximumBin()));
977  maxval = std::max(maxval, hdn->GetBinContent(hdn->GetMaximumBin()));
978  maxval = std::max(maxval, hsx->GetBinContent(hsx->GetMaximumBin()));
979 
980  double minval = nom->GetBinContent(nom->GetMinimumBin());
981  minval = std::min(minval, hup->GetBinContent(hup->GetMinimumBin()));
982  minval = std::min(minval, hdn->GetBinContent(hdn->GetMinimumBin()));
983  minval = std::min(minval, hsx->GetBinContent(hsx->GetMinimumBin()));
984 
985  if(maxval < minval)
986  std::swap(maxval, minval);
987 
988  double maxvalRat = hRat->GetBinContent(hRat->GetMaximumBin());
989  maxvalRat = std::max(maxvalRat, hUpRat->GetBinContent(hUpRat->GetMaximumBin()));
990  maxvalRat = std::max(maxvalRat, hDnRat->GetBinContent(hDnRat->GetMaximumBin()));
991  maxvalRat = std::max(maxvalRat, hSXRat->GetBinContent(hSXRat->GetMaximumBin()));
992 
993  double minvalRat = hRat->GetBinContent(hRat->GetMinimumBin());
994  for(int i = 1; i <= nBins; ++i) {
995  if (hUpRat->GetBinContent(i) > 0.)
996  minvalRat = std::min(minvalRat, hUpRat->GetBinContent(i));
997  if (hDnRat->GetBinContent(i) > 0.)
998  minvalRat = std::min(minvalRat, hDnRat->GetBinContent(i));
999  if (hSXRat->GetBinContent(i) > 0.)
1000  minvalRat = std::min(minvalRat, hSXRat->GetBinContent(i));
1001  }
1002 
1003  if (std::abs(maxvalRat - 1.) > std::abs(1. - minvalRat)) {
1004  minvalRat = 1 - std::abs(maxvalRat - 1.);
1005  maxvalRat = 1 + std::abs(maxvalRat - 1.);
1006  } else {
1007  minvalRat = 1 - std::abs(1. - minvalRat);
1008  maxvalRat = 1 + std::abs(1. - minvalRat);
1009  }
1010 
1011  CenterTitles(nom);
1012  nom->SetLineColor(kBlack);
1013  nom->SetMaximum(1.1*maxval);
1014  nom->SetMinimum(0);
1015  nom->SetTitle(fullTitle.c_str());
1016 
1017  CenterTitles(hup);
1018  hup->SetLineColor(kRed);
1019  hup->SetLineStyle(2);
1020  hup->SetMaximum(1.1*maxval);
1021  hup->SetMinimum(0);
1022  hup->SetTitle(fullTitle.c_str());
1023 
1024  CenterTitles(hdn);
1025  hdn->SetLineColor(kGreen+2);
1026  hdn->SetLineStyle(5);
1027  hdn->SetMaximum(1.1*maxval);
1028  hdn->SetMinimum(0);
1029  hdn->SetTitle(fullTitle.c_str());
1030 
1031  CenterTitles(hsx);
1032  hsx->SetLineColor(kBlue);
1033  hsx->SetLineStyle(2);
1034  hsx->SetMaximum(1.1*maxval);
1035  hsx->SetMinimum(0);
1036  hsx->SetTitle(fullTitle.c_str());
1037 
1038  CenterTitles(hRat);
1039  hRat->SetLineColor(kBlack);
1040  hRat->SetMaximum(1.1*maxvalRat);
1041  hRat->SetMinimum(0.9*minvalRat);
1042  hRat->SetTitle(fullTitleRat.c_str());
1043  hRat->GetYaxis()->SetTitleSize(.06);
1044  hRat->GetXaxis()->SetTitleSize(.06);
1045  hRat->GetXaxis()->SetLabelSize(.06);
1046  hRat->GetYaxis()->SetLabelSize(.06);
1047  hRat->GetYaxis()->SetTitleOffset(.5);
1048  hRat->GetXaxis()->SetTitleOffset(.8);
1049 
1050  CenterTitles(hUpRat);
1051  hUpRat->SetLineColor(kRed);
1052  hUpRat->SetLineStyle(2);
1053  hUpRat->SetMaximum(1.1*maxvalRat);
1054  hUpRat->SetMinimum(0.9*minvalRat);
1055  hUpRat->SetTitle(fullTitleRat.c_str());
1056  hUpRat->GetYaxis()->SetTitleSize(.06);
1057  hUpRat->GetXaxis()->SetTitleSize(.06);
1058  hUpRat->GetXaxis()->SetLabelSize(.06);
1059  hUpRat->GetYaxis()->SetLabelSize(.06);
1060  hUpRat->GetYaxis()->SetTitleOffset(.5);
1061  hUpRat->GetXaxis()->SetTitleOffset(.8);
1062 
1063  CenterTitles(hDnRat);
1064  hDnRat->SetLineColor(kGreen+2);
1065  hDnRat->SetLineStyle(5);
1066  hDnRat->SetMaximum(1.1*maxvalRat);
1067  hDnRat->SetMinimum(0.9*minvalRat);
1068  hDnRat->SetTitle(fullTitleRat.c_str());
1069  hDnRat->GetYaxis()->SetTitleSize(.06);
1070  hDnRat->GetXaxis()->SetTitleSize(.06);
1071  hDnRat->GetXaxis()->SetLabelSize(.06);
1072  hDnRat->GetYaxis()->SetLabelSize(.06);
1073  hDnRat->GetYaxis()->SetTitleOffset(.5);
1074  hDnRat->GetXaxis()->SetTitleOffset(.8);
1075 
1076  CenterTitles(hSXRat);
1077  hSXRat->SetLineColor(kBlue);
1078  hSXRat->SetLineStyle(2);
1079  hSXRat->SetMaximum(1.1*maxvalRat);
1080  hSXRat->SetMinimum(0.9*minvalRat);
1081  hSXRat->SetTitle(fullTitleRat.c_str());
1082  hSXRat->GetYaxis()->SetTitleSize(.09);
1083  hSXRat->GetXaxis()->SetTitleSize(.06);
1084  hSXRat->GetXaxis()->SetLabelSize(.06);
1085  hSXRat->GetYaxis()->SetLabelSize(.06);
1086  hSXRat->GetYaxis()->SetTitleOffset(.5);
1087  hSXRat->GetXaxis()->SetTitleOffset(.8);
1088 
1089  double xL = 0.6, xR = 0.85;
1090  double yB = 0.6, yT = 0.85;
1091  TLegend* leg = new TLegend(xL, yB, xR, yT);
1092  leg->AddEntry(nom, "Nominal", "l");
1093  leg->AddEntry(hup, "Flat Scale Up", "l");
1094  leg->AddEntry(hdn, "Flat Scale Down", "l");
1095  leg->AddEntry(hsx, "Calibration shape", "l");
1096 
1097  TCanvas* c = new TCanvas(name.c_str(), title.c_str(), 800, 800);
1098  TPad* pSpecs = new TPad("pSpecs", "", 0., 0.375, 1., 1.);
1099  TPad* pRatio = new TPad("pRatio", "", 0., 0., 1., 0.375);
1100  pSpecs->Draw();
1101  pRatio->Draw();
1102 
1103  pSpecs->cd();
1104  nom->Draw("hist");
1105  hup->Draw("hist same");
1106  hdn->Draw("hist same");
1107  hsx->Draw("hist same");
1108  leg->Draw();
1109  Simulation();
1110 
1111  pRatio->cd();
1112  hRat->Draw("hist");
1113  hUpRat->Draw("hist same");
1114  hDnRat->Draw("hist same");
1115  hSXRat->Draw("hist same");
1116 
1117  c->Update();
1118  rootOUT->WriteTObject(c);
1119 
1120  return;
1121 
1122  }
1123 
1124  //----------------------------------------------------------------------------
1125  /// Plot multiple systematics
1126  void PlotMultiSyst(TH1* h, std::vector<TH1*> hp,
1127  TDirectory* out, FILE* text, strings strs) {
1128 
1129  std::string xlabel = strs.fXLabel;
1130  std::string ylabel = "Events / " + strs.fPOT + " / 0.25 GeV";
1131  std::string title = strs.fDet + " " + strs.fComponent +
1132  " Error for " + strs.fSystL + " Systematic";
1133  std::string fullTitle = ";" + xlabel + ";" + ylabel;
1134  std::string name = "c" + strs.fComponent + strs.fDet +
1135  strs.fSystS + strs.fSample;
1136  std::string ylabelRat = "Shifted / Nominal";
1137  std::string fullTitleRat = ";" + xlabel + ";" + ylabelRat;
1138 
1139  //
1140  TH1* hCVF = hp.at(0);
1141  TH1* hCV = (TH1*)h->Clone();
1142  std::map<int, float> cv_events_in_bins;
1143  std::map<int, std::vector<float> > univ_events_in_bins;
1144  unsigned int count = 0;
1145 
1146  for (auto& it: hp) {
1147  if (count!= hp.size()) {
1148  TH1* curhist = it;
1149  int nbinsx = curhist->GetNbinsX();
1150  for (int i = 1; i <= nbinsx; i++)
1151  univ_events_in_bins[i].push_back(curhist->GetBinContent(i));
1152  }
1153  else if (count== hp.size()) {
1154  int nbinsx = hCV->GetNbinsX();
1155  for (int i = 1; i <= nbinsx; i++)
1156  cv_events_in_bins[i] = hCV->GetBinContent(i);
1157  }
1158  count++;
1159  }
1160 
1161  TH1* hUp = (TH1*)hCVF->Clone("hUp");
1162  TH1* hLow = (TH1*)hCVF->Clone("hLow");
1163  int binIt = 1;
1164 
1165  for (auto& it: univ_events_in_bins) {
1166  double cv_val = hCVF->GetBinContent( it.first);
1167  double err = 0.0;
1168  for (unsigned int iuniv=0;iuniv<hp.size();iuniv++) {
1169  err += pow(it.second[iuniv] - cv_val,2);
1170  }
1171  err /= double(hp.size());
1172  err = sqrt(err);
1173  hUp->SetBinContent(binIt, cv_val+err);
1174  hLow->SetBinContent(binIt,cv_val-err);
1175  binIt++;
1176  }
1177  TH1* hRat = (TH1*)h->Clone();
1178  TH1* hUpRat=(TH1*)hUp->Clone();
1179  TH1* hLowRat=(TH1*)hLow->Clone();
1180  int nBins = h->GetNbinsX();
1181  double percent_diff = 100. * std::abs(hUp->Integral(1, nBins) - h->Integral(1, nBins))/
1182  h->Integral(1, nBins);
1183  double percent_diff1 = 100. * std::abs(hLow->Integral(1, nBins) - h->Integral(1, nBins))/
1184  h->Integral(1, nBins);
1185 
1186  fprintf(text, "Up %s %s %.2s: +/-%6.4f%%\n",
1187  strs.fSystS.c_str(), strs.fComponent.c_str(), strs.fDet.c_str(),
1188  percent_diff);
1189  fprintf(text, "Low %s %s %.2s: +/-%6.4f%%\n",
1190  strs.fSystS.c_str(), strs.fComponent.c_str(), strs.fDet.c_str(),
1191  percent_diff1);
1192 
1193  for (int i = 0; i <= nBins+1; ++i) {
1194  hRat->SetBinContent(i, 1.);
1195  }
1196  hUpRat->Divide(h);
1197  hLowRat->Divide(h);
1198 
1199  double maxval = h->GetBinContent(h->GetMaximumBin());
1200  maxval = std::max(maxval, hUp->GetBinContent(hUp->GetMaximumBin()));
1201 
1202  double minval = h->GetBinContent(h->GetMinimumBin());
1203  minval = std::min(minval, hLow->GetBinContent(hLow->GetMinimumBin()));
1204 
1205  if(maxval < minval) { std::swap(maxval, minval); }
1206 
1207  double maxvalRat = hRat->GetBinContent(hRat->GetMaximumBin());
1208  maxvalRat = std::max(maxvalRat, hUpRat->GetBinContent(hUpRat->GetMaximumBin()));
1209 
1210  double minvalRat = hLowRat->GetBinContent(hLowRat->GetMinimumBin());
1211  for (int i = 1; i <= nBins; ++i)
1212  if(hUpRat->GetBinContent(i) > 0.)
1213  minvalRat = std::min(minvalRat, hLowRat->GetBinContent(i));
1214 
1215  if (std::abs(maxvalRat - 1.) > std::abs(1. - minvalRat)) {
1216  minvalRat = 1 - std::abs(maxvalRat - 1.);
1217  maxvalRat = 1 + std::abs(maxvalRat - 1.);
1218  } else {
1219  minvalRat = 1 - std::abs(1. - minvalRat);
1220  maxvalRat = 1 + std::abs(1. - minvalRat);
1221  }
1222 
1223  CenterTitles(h);
1224  h->SetLineColor(kBlack);
1225  h->SetMaximum(1.1*maxval);
1226  h->SetMinimum(0);
1227  h->SetTitle(fullTitle.c_str());
1228 
1229  CenterTitles(hLow);
1230  hLow->SetLineColor(kRed);
1231  hLow->SetLineStyle(2);
1232  hLow->SetMaximum(1.1*maxval);
1233  hLow->SetMinimum(0);
1234  hLow->SetTitle(fullTitle.c_str());
1235 
1236  CenterTitles(hUp);
1237  hUp->SetLineColor(kGreen+2);
1238  hUp->SetLineStyle(2);
1239  hUp->SetMaximum(1.1*maxval);
1240  hUp->SetMinimum(0);
1241  hUp->SetTitle(fullTitle.c_str());
1242 
1243  CenterTitles(hRat);
1244  hRat->SetLineColor(kBlack);
1245  hRat->SetMaximum(1.1*maxvalRat);
1246  hRat->SetMinimum(0.9*minvalRat);
1247  hRat->SetTitle(fullTitleRat.c_str());
1248 
1249  CenterTitles(hUpRat);
1250  hUpRat->SetLineColor(kGreen+2);
1251  hUpRat->SetLineStyle(2);
1252  hUpRat->SetMaximum(1.1*maxvalRat);
1253  hUpRat->SetMinimum(0.9*minvalRat);
1254 
1255  CenterTitles(hLowRat);
1256  hLowRat->SetLineColor(kRed);
1257  hLowRat->SetLineStyle(2);
1258  hLowRat->SetMaximum(1.1*maxvalRat);
1259  hLowRat->SetMinimum(0.9*minvalRat);
1260  double xL = 0.6, xR = 0.85;
1261  double yB = 0.6, yT = 0.85;
1262  TLegend* leg = new TLegend(xL, yB, xR, yT);
1263  leg->AddEntry(h, "Nominal", "l");
1264  leg->AddEntry(hUp, "Shifted Up 34%", "l");
1265  leg->AddEntry(hLow, "Shifted down 34%", "l");
1266 
1267  TCanvas* c = new TCanvas(name.c_str(), title.c_str(), 800, 800);
1268  TPad* pSpecs = new TPad("pSpecs", "", 0., 0.375, 1., 1.);
1269  TPad* pRatio = new TPad("pRatio", "", 0., 0., 1., 0.375);
1270  pSpecs->Draw();
1271  pRatio->Draw();
1272 
1273  pSpecs->cd();
1274  h->Draw("hist");
1275  hUp->Draw("hist same");
1276  hLow->Draw("hist same");
1277  c->Update();
1278  leg->Draw();
1279  Simulation();
1280 
1281  pRatio->cd();
1282  hRat->Draw("hist");
1283  hUpRat->Draw("hist same");
1284  hLowRat->Draw("hist same");
1285  c->Update();
1286  out->WriteTObject(c);
1287 
1288  return;
1289 
1290  }
1291 
1292  //----------------------------------------------------------------------------
1293  /// Plot multiple systematics
1294  void PlotMultiSyst(TH1* nom, std::vector<NDPredictionSterile*> shifts,
1295  TDirectory* out, FILE* text, strings strs,
1296  double POT, bool NC, bool use_max) {
1297 
1298  std::vector<TH1*> hp;
1299  for(auto shift: shifts)
1300  hp.push_back(NC ? GetNC(shift, POT) : GetBG(shift, POT));
1301 
1302  PlotMultiSyst(nom, hp, out, text, strs);
1303  return;
1304 
1305  }
1306 
1307  //----------------------------------------------------------------------------
1308  /// Plot multiple systematics
1309  void PlotMultiSyst(TH1* nom, std::vector<FDPredictionSterile*> shifts,
1310  TDirectory* out, FILE* text, strings strs,
1311  double POT, bool NC, bool use_max) {
1312 
1313  std::vector<TH1*> hp;
1314  for(auto shift: shifts)
1315  hp.push_back(NC ? GetNC(shift, POT) : GetBG(shift, POT));
1316 
1317  PlotMultiSyst(nom, hp, out, text, strs);
1318  return;
1319 
1320  }
1321 
1322  //----------------------------------------------------------------------------
1323  /// Plot multiple systematics
1324  void PlotMultiSyst(TH1* nom, std::vector<PredictionSterile*> shifts,
1325  TDirectory* out, FILE* text, strings strs,
1326  double POT, bool NC, bool use_max) {
1327 
1328  std::vector<TH1*> hp;
1329  for(auto shift: shifts)
1330  hp.push_back(NC ? GetNC(shift, POT) : GetBG(shift, POT));
1331 
1332  PlotMultiSyst(nom, hp, out, text, strs);
1333  return;
1334 
1335  }
1336 
1337  //----------------------------------------------------------------------------
1338  /// Plot FD syst shifted and nominal spectra
1339  /// This is for plotting an on/off style shift
1340  void PlotSyst(TH1* nom, IPrediction* shifts,
1341  TDirectory* out, FILE* text, strings strs,
1342  double POT, bool NC, bool use_max) {
1343 
1344  TH1* hp1 = (NC ? GetNC(shifts, POT) : GetBG(shifts, POT));
1345 
1346  AddErrorInQuadrature(nom, hp1, use_max);
1347  PlotSyst(nom, hp1, out, text, strs);
1348 
1349  return;
1350 
1351  }
1352 
1353  //----------------------------------------------------------------------------
1354  /// Plot FD syst shifted and nominal spectra
1355  /// This is for plotting a syst with multiple sigma levels
1356  void PlotSyst(TH1* nom, std::map<int, IPrediction*> shifts,
1357  TDirectory* out, FILE* text, strings strs,
1358  double POT, bool NC, bool use_max) {
1359 
1360  TH1* hp1 = (NC ? GetNC(shifts[+1], POT) : GetBG(shifts[+1], POT));
1361  TH1* hm1 = (NC ? GetNC(shifts[-1], POT) : GetBG(shifts[-1], POT));
1362  TH1 *hp2, *hm2;
1363  if (shifts.count(+2))
1364  hp2 = (NC ? GetNC(shifts[+2], POT) : GetBG(shifts[+2], POT));
1365  if (shifts.count(-2))
1366  hm2 = (NC ? GetNC(shifts[-2], POT) : GetBG(shifts[-2], POT));
1367 
1368  AddErrorInQuadrature(nom, hp1, hm1, use_max);
1369  if (hp1 and hm1 and hp2 and hm2)
1370  PlotSyst(nom, hp1, hm1, hp2, hm2, out, text, strs);
1371  else if (hp1 and hm1)
1372  PlotSyst(nom, hp1, hm1, out, text, strs);
1373 
1374  return;
1375 
1376  }
1377 
1378  //----------------------------------------------------------------------------
1379  void PlotSyst(TH1* nom, std::map<int, NDPredictionSterile*> shifts,
1380  TDirectory* out, FILE* text, strings strs,
1381  double POT, bool NC, bool use_max) {
1382 
1383  TH1* hp1 = (NC ? GetNC(shifts[+1], POT) : GetBG(shifts[+1], POT));
1384  TH1* hm1 = (NC ? GetNC(shifts[-1], POT) : GetBG(shifts[-1], POT));
1385  TH1 *hp2, *hm2;
1386  if (shifts.count(+2))
1387  hp2 = (NC ? GetNC(shifts[+2], POT) : GetBG(shifts[+2], POT));
1388  if (shifts.count(-2))
1389  hm2 = (NC ? GetNC(shifts[-2], POT) : GetBG(shifts[-2], POT));
1390 
1391  AddErrorInQuadrature(nom, hp1, hm1, use_max);
1392  if (hp1 and hm1 and hp2 and hm2)
1393  PlotSyst(nom, hp1, hm1, hp2, hm2, out, text, strs);
1394  else if (hp1 and hm1)
1395  PlotSyst(nom, hp1, hm1, out, text, strs);
1396 
1397  return;
1398 
1399  }
1400 
1401  //----------------------------------------------------------------------------
1402  void PlotSyst(TH1* nom, std::map<int, FDPredictionSterile*> shifts,
1403  TDirectory* out, FILE* text, strings strs,
1404  double POT, bool NC, bool use_max) {
1405 
1406  TH1* hp1 = (NC ? GetNC(shifts[+1], POT) : GetBG(shifts[+1], POT));
1407  TH1* hm1 = (NC ? GetNC(shifts[-1], POT) : GetBG(shifts[-1], POT));
1408  TH1 *hp2, *hm2;
1409  if (shifts.count(+2))
1410  hp2 = (NC ? GetNC(shifts[+2], POT) : GetBG(shifts[+2], POT));
1411  if (shifts.count(-2))
1412  hm2 = (NC ? GetNC(shifts[-2], POT) : GetBG(shifts[-2], POT));
1413 
1414  AddErrorInQuadrature(nom, hp1, hm1, use_max);
1415  if (hp1 and hm1 and hp2 and hm2)
1416  PlotSyst(nom, hp1, hm1, hp2, hm2, out, text, strs);
1417  else if (hp1 and hm1)
1418  PlotSyst(nom, hp1, hm1, out, text, strs);
1419 
1420  return;
1421 
1422  }
1423 
1424  //----------------------------------------------------------------------------
1425  void PlotSyst(TH1* nom, std::map<int, PredictionSterile*> shifts,
1426  TDirectory* out, FILE* text, strings strs,
1427  double POT, bool NC, bool use_max) {
1428 
1429  TH1* hp1 = (NC ? GetNC(shifts[+1], POT) : GetBG(shifts[+1], POT));
1430  TH1* hm1 = (NC ? GetNC(shifts[-1], POT) : GetBG(shifts[-1], POT));
1431  TH1 *hp2, *hm2;
1432  if (shifts.count(+2))
1433  hp2 = (NC ? GetNC(shifts[+2], POT) : GetBG(shifts[+2], POT));
1434  if (shifts.count(-2))
1435  hm2 = (NC ? GetNC(shifts[-2], POT) : GetBG(shifts[-2], POT));
1436 
1437  AddErrorInQuadrature(nom, hp1, hm1, use_max);
1438  if (hp1 and hm1 and hp2 and hm2)
1439  PlotSyst(nom, hp1, hm1, hp2, hm2, out, text, strs);
1440  else if (hp1 and hm1)
1441  PlotSyst(nom, hp1, hm1, out, text, strs);
1442 
1443  return;
1444 
1445  }
1446 
1447  //----------------------------------------------------------------------------
1448  /// Plot the systematic error band around a nominal spectrum
1449  /// The errors of the input histogram should be squared
1450  void PlotSystBand(TH1* h, TDirectory* out, FILE* text, strings strs)
1451  {
1452  // Fix the errors by taking the square root!
1453  SqrtError(h);
1454 
1455  std::string xlabel = strs.fXLabel;
1456  std::string ylabel = "Events / " + strs.fPOT + " / 0.25 GeV";
1457  std::string title = strs.fDet + " " + strs.fComponent +
1458  " Error Band for All " + strs.fSystType + " Systematics";
1459  std::string fullTitle = ";" + xlabel + ";" + ylabel;
1460  std::string name = "c" + strs.fComponent + strs.fDet +
1461  strs.fSystType + "Systs" + strs.fSample;
1462  std::string ylabelRat = "Shifted / Nominal";
1463  std::string fullTitleRat = ";" + xlabel + ";" + ylabelRat;
1464 
1465  int nBins = h->GetNbinsX();
1466  double maxval = 0.;
1467  double nom_evts = h->Integral(1, nBins);
1468  double evts_diff = 0.;
1469  double percent_diff = 100.;
1470  for(int i = 1; i <= nBins; ++i) {
1471  maxval = std::max(maxval, h->GetBinContent(i) + h->GetBinError(i));
1472  evts_diff += std::abs(h->GetBinError(i));
1473  }
1474  percent_diff *= (evts_diff/nom_evts);
1475 
1476  fprintf(text, "\nAll systs %.2s %.2s: +/-%6.4f%%\n",
1477  strs.fComponent.c_str(), strs.fDet.c_str(),
1478  percent_diff);
1479 
1480  CenterTitles(h);
1481  h->SetLineColor(kTotalMCColor);
1482  h->SetMaximum(1.1*maxval);
1483  h->SetMinimum(0);
1484  h->SetTitle(fullTitle.c_str());
1485 
1486  TH1* hErr = (TH1*)h->Clone();
1487  CenterTitles(hErr);
1488  hErr->SetFillColor(kTotalMCColor-9);
1489  hErr->SetMaximum(1.1*maxval);
1490  hErr->SetMinimum(0);
1491 
1492  double maxvalRat = 1.;
1493  TH1* hRat = (TH1*)h->Clone();
1494  TH1* hUpp = (TH1*)h->Clone();
1495  TH1* hDwn = (TH1*)h->Clone();
1496  for(int i = 1, n = hRat->GetNbinsX(); i <= n; ++i) {
1497  double binVal = hRat->GetBinContent(i);
1498  double errVal = hRat->GetBinError(i);
1499  if(binVal <= 0.) {
1500  binVal = 1.;
1501  errVal = 0.;
1502  }
1503 
1504  hRat->SetBinContent(i, 1);
1505  hUpp->SetBinContent(i, (binVal + errVal)/binVal);
1506  hDwn->SetBinContent(i, (binVal - errVal)/binVal);
1507 
1508  double diff = (binVal + errVal)/binVal;
1509  if(diff > maxvalRat) {
1510  maxvalRat = diff;
1511  }
1512  }
1513  double minvalRat = 1. - (maxvalRat-1.);
1514 
1515  CenterTitles(hRat);
1516  hRat->SetLineColor(kBlack);
1517  hRat->SetMaximum(1.05*maxvalRat);
1518  hRat->SetMinimum(0.95*minvalRat);
1519  hRat->SetTitle(fullTitleRat.c_str());
1520 
1521  CenterTitles(hUpp);
1522  hUpp->SetLineColor(kRed);
1523  hUpp->SetLineStyle(2);
1524  hUpp->SetMaximum(1.05*maxvalRat);
1525  hUpp->SetMinimum(0.95*minvalRat);
1526  hUpp->SetTitle(fullTitleRat.c_str());
1527 
1528  CenterTitles(hDwn);
1529  hDwn->SetLineColor(kRed);
1530  hDwn->SetLineStyle(2);
1531  hDwn->SetMaximum(1.05*maxvalRat);
1532  hDwn->SetMinimum(0.95*minvalRat);
1533  hDwn->SetTitle(fullTitleRat.c_str());
1534 
1535  TCanvas* c = new TCanvas(name.c_str(), title.c_str(), 800, 800);
1536  TPad* pSpecs = new TPad("pSpecs", "", 0., 0.375, 1., 1.);
1537  TPad* pRatio = new TPad("pRatio", "", 0., 0., 1., 0.375);
1538  pSpecs->Draw();
1539  pRatio->Draw();
1540 
1541  pSpecs->cd();
1542  hErr->Draw("e2");
1543  h->Draw("hist same");
1544  Simulation();
1545 
1546  pRatio->cd();
1547  hRat->Draw("hist");
1548  hUpp->Draw("hist same");
1549  hDwn->Draw("hist same");
1550 
1551  c->Update();
1552  out->WriteTObject(c);
1553  return;
1554  }
1555 
1556  //----------------------------------------------------------------------------
1557  void SaveMaps(TDirectory* out,
1558  std::map<std::string, NDPredictionSterile*> predsND_nominal,
1559  std::map<std::string, std::map<std::string, std::map<int, NDPredictionSterile*> > > predsND_shifted) {
1560 
1561  std::string dir, sep = "__";
1562 
1563  for (const auto& pred : predsND_nominal) {
1564  dir = "NDPred" + sep + pred.first;
1565  pred.second->SaveTo(out->mkdir(dir.c_str()));
1566  }
1567 
1568  for (const auto& sample : predsND_shifted)
1569  for (const auto& shiftlabel : sample.second)
1570  for (const auto& pred : shiftlabel.second) {
1571  dir = Form("NDPred%s%s%s%s%s%d", sep.c_str(),
1572  sample.first.c_str(), sep.c_str(),
1573  shiftlabel.first.c_str(), sep.c_str(), pred.first);
1574  pred.second->SaveTo(out->mkdir(dir.c_str()));
1575  }
1576 
1577  return;
1578 
1579  }
1580 
1581  //----------------------------------------------------------------------------
1582  void SaveMaps(TDirectory* out,
1583  std::map<std::string, FDPredictionSterile*> predsFD_nominal,
1584  std::map<std::string, std::map<std::string, std::map<int, FDPredictionSterile*> > > predsFD_shifted) {
1585 
1586  std::string dir, sep = "__";
1587 
1588  for (const auto& pred : predsFD_nominal) {
1589  dir = "FDPred" + sep + pred.first;
1590  pred.second->SaveTo(out->mkdir(dir.c_str()));
1591  }
1592 
1593  for (const auto& sample : predsFD_shifted)
1594  for (const auto& shiftlabel : sample.second)
1595  for (const auto& pred : shiftlabel.second) {
1596  dir = Form("FDPred%s%s%s%s%s%d", sep.c_str(),
1597  sample.first.c_str(), sep.c_str(),
1598  shiftlabel.first.c_str(), sep.c_str(), pred.first);
1599  pred.second->SaveTo(out->mkdir(dir.c_str()));
1600  }
1601 
1602  return;
1603 
1604  }
1605 
1606  //----------------------------------------------------------------------------
1607  void SaveMaps(TDirectory* out,
1608  std::map<std::string, PredictionSterile*> preds_nominal,
1609  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > > preds_shifted) {
1610 
1611  std::string dir, sep = "__";
1612 
1613  for (const auto& pred : preds_nominal) {
1614  dir = "ExtrapPred" + sep + pred.first;
1615  pred.second->SaveTo(out->mkdir(dir.c_str()));
1616  }
1617 
1618  for (const auto& sample : preds_shifted)
1619  for (const auto& shiftlabel : sample.second)
1620  for (const auto& pred : shiftlabel.second) {
1621  dir = Form("ExtrapPred%s%s%s%s%s%d", sep.c_str(),
1622  sample.first.c_str(), sep.c_str(),
1623  shiftlabel.first.c_str(), sep.c_str(), pred.first);
1624  pred.second->SaveTo(out->mkdir(dir.c_str()));
1625  }
1626 
1627  return;
1628 
1629  }
1630 
1631 }
static std::unique_ptr< FDPredictionSterile > LoadFrom(TDirectory *dir)
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
osc::OscCalculatorDumb calc
Oscillation analysis framework, runs over CAF files outside of ART.
void PlotMultiSyst(TH1 *h, std::vector< TH1 * > hp, TDirectory *out, FILE *text, strings strs)
Plot multiple systematics.
Definition: PPFXHelper.h:464
set< int >::iterator it
std::string fComponent
Definition: PPFXHelper.h:103
Adapt the PMNS_Sterile calculator to standard interface.
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:553
virtual Spectrum PredictComponent(osc::IOscCalculator *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
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
constexpr T pow(T x)
Definition: pow.h:75
static std::unique_ptr< NDPredictionSterile > LoadFrom(TDirectory *dir)
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
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:1128
void Simulation()
Definition: nue_pid_effs.C:49
float abs(float number)
Definition: d0nt_math.hpp:39
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:37
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
base_types push_back(int_type())
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
static std::unique_ptr< PredictionSterile > LoadFrom(TDirectory *dir)
std::string fPOT
Definition: PPFXHelper.h:105
float bin[41]
Definition: plottest35.C:14
const char sep
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
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
Int_t nbinsx
Definition: plot.C:23
TDirectory * dir
Definition: macro.C:5
std::string fSystType
Definition: PPFXHelper.h:107
Neutral-current interactions.
Definition: IPrediction.h:40
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
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:17
osc::OscCalculatorSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
Definition: Calcs.cxx:97
T min(const caf::Proxy< T > &a, T b)
All neutrinos, any flavor.
Definition: IPrediction.h:26
TH1 * GetBG(IDecomp *specs, double POT)
Definition: PPFXHelper.h:236
float count
Definition: make_cosmics.C:27
A helper structure to contain a group of string for plotting.
Definition: PPFXHelper.h:101
std::string fSample
Definition: PPFXHelper.h:106
h
Definition: demo3.py:41