PPFXHelper.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 
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"
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  /// \param F
464  void PlotMultiSyst(TH1* h, std::vector<TH1*> hp, TDirectory* out, FILE* text, strings strs)
465  {
466  // Construct histogram names and titles
467  std::string xlabel = strs.fXLabel;
468  std::string ylabel = "Events / " + strs.fPOT + " / 0.25 GeV";
469  std::string title = strs.fDet + " " + strs.fComponent +
470  " Error for " + strs.fSystL + " Systematic";
471  std::string fullTitle = ";" + xlabel + ";" + ylabel;
472 // std::string fullTitle = title + ";" + xlabel + ";" + ylabel;
473  std::string name = "c" + strs.fComponent + strs.fDet +
474  strs.fSystS + strs.fSample;
475  std::string ylabelRat = "Shifted / Nominal";
476  std::string fullTitleRat = ";" + xlabel + ";" + ylabelRat;
477  /* TH1* hRat;
478  std::vector<TH1*>hp1V;
479  std::vector<TH1*>hp1RatV;
480  hp1V.clear();
481  hp1RatV.clear();
482  */ // for(const auto hp1 :hp){
483  TH1* hCVF = hp.at(0);
484  TH1* hCV = (TH1*)h->Clone();
485  map< int,float > cv_events_in_bins;
486  map<int, vector<float> > univ_events_in_bins;
487  unsigned int count = 0;
488 
489  for(auto& it: hp){
490  if(count!= hp.size()){
491  TH1* curhist = it;
492  int nbinsx = curhist->GetNbinsX();
493  for(int i = 1; i <= nbinsx; i++){
494  univ_events_in_bins[i].push_back(curhist->GetBinContent(i));
495  }
496  }
497  else if(count== hp.size()){
498  int nbinsx = hCV->GetNbinsX();
499  for(int i = 1; i <= nbinsx; i++){
500  cv_events_in_bins[i] = hCV->GetBinContent(i);
501  }
502  }
503  count++;
504  }
505 
506  TH1* hUp = (TH1*)hCVF->Clone("hUp");
507  TH1* hLow = (TH1*)hCVF->Clone("hLow");
508  int binIt = 1;
509  for(auto& it: univ_events_in_bins){
510  double cv_val = hCVF->GetBinContent( it.first);
511  double err = 0.0;
512  for(unsigned int iuniv=0;iuniv<hp.size();iuniv++){
513  err += pow(it.second[iuniv] - cv_val,2);
514  }
515  err /= double(hp.size());
516  err = sqrt(err);
517  hUp->SetBinContent(binIt, cv_val+err);
518  hLow->SetBinContent(binIt,cv_val-err);
519  binIt++;
520  }
521 
522 // TH1* hp1=hp.at(0);
523  TH1* hRat = (TH1*)h->Clone();
524  // TH1* hp1Rat = (TH1*)hp1->Clone();
525  TH1* hUpRat=(TH1*)hUp->Clone();
526  TH1* hLowRat=(TH1*)hLow->Clone();
527  int nBins = h->GetNbinsX();
528  double percent_diff = 100. * std::abs(hUp->Integral(1, nBins) - h->Integral(1, nBins))/
529  h->Integral(1, nBins);
530  double percent_diff1 = 100. * std::abs(hLow->Integral(1, nBins) - h->Integral(1, nBins))/
531  h->Integral(1, nBins);
532 
533  fprintf(text, "Up %s %s %.2s: +/-%6.4f%%\n",
534  strs.fSystS.c_str(), strs.fComponent.c_str(), strs.fDet.c_str(),
535  percent_diff);
536  fprintf(text, "Low %s %s %.2s: +/-%6.4f%%\n",
537  strs.fSystS.c_str(), strs.fComponent.c_str(), strs.fDet.c_str(),
538  percent_diff1);
539 
540 
541  for(int i = 0; i <= nBins+1; ++i) {
542  hRat->SetBinContent(i, 1.);
543  }
544  hUpRat->Divide(h);
545  hLowRat->Divide(h);
546 
547  double maxval = h->GetBinContent(h->GetMaximumBin());
548  maxval = std::max(maxval, hUp->GetBinContent(hUp->GetMaximumBin()));
549 
550  double minval = h->GetBinContent(h->GetMinimumBin());
551  minval = std::min(minval, hLow->GetBinContent(hLow->GetMinimumBin()));
552 
553  if(maxval < minval) { std::swap(maxval, minval); }
554 
555  double maxvalRat = hRat->GetBinContent(hRat->GetMaximumBin());
556  maxvalRat = std::max(maxvalRat, hUpRat->GetBinContent(hUpRat->GetMaximumBin()));
557 
558  double minvalRat = hLowRat->GetBinContent(hLowRat->GetMinimumBin());
559  for(int i = 1; i <= nBins; ++i) {
560  if(hUpRat->GetBinContent(i) > 0.) {
561  minvalRat = std::min(minvalRat, hLowRat->GetBinContent(i));
562  }
563  }
564 
565  if(std::abs(maxvalRat - 1.) > std::abs(1. - minvalRat)) {
566  minvalRat = 1 - std::abs(maxvalRat - 1.);
567  maxvalRat = 1 + std::abs(maxvalRat - 1.);
568  }
569  else {
570  minvalRat = 1 - std::abs(1. - minvalRat);
571  maxvalRat = 1 + std::abs(1. - minvalRat);
572  }
573 
574  CenterTitles(h);
575  h->SetLineColor(kBlack);
576  h->SetMaximum(1.1*maxval);
577  h->SetMinimum(0);
578  h->SetTitle(fullTitle.c_str());
579 
580  CenterTitles(hLow);
581  hLow->SetLineColor(kRed);
582  hLow->SetLineStyle(2);
583  hLow->SetMaximum(1.1*maxval);
584  hLow->SetMinimum(0);
585  hLow->SetTitle(fullTitle.c_str());
586 
587  CenterTitles(hUp);
588  hUp->SetLineColor(kGreen+2);
589  hUp->SetLineStyle(2);
590  hUp->SetMaximum(1.1*maxval);
591  hUp->SetMinimum(0);
592  hUp->SetTitle(fullTitle.c_str());
593 
594 
595  CenterTitles(hRat);
596  hRat->SetLineColor(kBlack);
597  hRat->SetMaximum(1.1*maxvalRat);
598  hRat->SetMinimum(0.9*minvalRat);
599  hRat->SetTitle(fullTitleRat.c_str());
600 
601  CenterTitles(hUpRat);
602  hUpRat->SetLineColor(kGreen+2);
603  hUpRat->SetLineStyle(2);
604  hUpRat->SetMaximum(1.1*maxvalRat);
605  hUpRat->SetMinimum(0.9*minvalRat);
606 
607  CenterTitles(hLowRat);
608  hLowRat->SetLineColor(kRed);
609  hLowRat->SetLineStyle(2);
610  hLowRat->SetMaximum(1.1*maxvalRat);
611  hLowRat->SetMinimum(0.9*minvalRat);
612  //hp1Rat->SetTitle(fullTitleRat.c_str());
613  // hp1RatV.push_back(hp1Rat);
614  // hp1V.push_back(hp1);
615 //}
616  double xL = 0.6, xR = 0.85;
617  double yB = 0.6, yT = 0.85;
618  TLegend* leg = new TLegend(xL, yB, xR, yT);
619  leg->AddEntry(h, "Nominal", "l");
620  leg->AddEntry(hUp, "Shifted Up 34%", "l");
621  leg->AddEntry(hLow, "Shifted down 34%", "l");
622 
623  // for(const auto hp1 :hp){
624  TCanvas* c = new TCanvas(name.c_str(), title.c_str(), 800, 800);
625  TPad* pSpecs = new TPad("pSpecs", "", 0., 0.375, 1., 1.);
626  TPad* pRatio = new TPad("pRatio", "", 0., 0., 1., 0.375);
627  pSpecs->Draw();
628  pRatio->Draw();
629 
630  pSpecs->cd();
631  h->Draw("hist");
632  hUp->Draw("hist same");
633  hLow->Draw("hist same");
634  c->Update();
635  leg->Draw();
636  Simulation();
637 
638  pRatio->cd();
639  hRat->Draw("hist");
640 // for(auto hp1Rat:hp1RatV){
641  hUpRat->Draw("hist same");
642  hLowRat->Draw("hist same");
643 // }
644  c->Update();
645 //}
646  out->WriteTObject(c);
647  // }
648  return;
649  }
650 
651  //----------------------------------------------------------------------------
652  /// Plot a nominal spectrum and a shifted spectra
653  /// This is for a systematic with multiple sigma levels
654  /// \param h The nominal spectrum
655  /// \param hp1 The shifted spectrum, +1 sigma
656  /// \param hm1 The shifted spectrum, -1 sigma
657  /// \param hp2 The shifted spectrum, +2 sigma
658  /// \param hm2 The shifted spectrum, -2 sigma
659  /// \param out The directory/file for root output
660  /// \param text The file for text output
661  /// \param strs Contains strings for pretty plots
662  void PlotSyst(TH1* h, TH1* hp1, TH1* hm1, TH1* hp2, TH1* hm2,
663  TDirectory* out, FILE* text, strings strs)
664  {
665  std::string xlabel = strs.fXLabel;
666  std::string ylabel = "Events / " + strs.fPOT + " / 0.25 GeV";
667  std::string title = strs.fDet + " " + strs.fComponent +
668  " Error for " + strs.fSystL + " Systematic";
669  std::string fullTitle = ";" + xlabel + ";" + ylabel;
670 // std::string fullTitle = title + ";" + xlabel + ";" + ylabel;
671  std::string name = "c" + strs.fComponent + strs.fDet +
672  strs.fSystS + strs.fSample;
673  std::string ylabelRat = "Shifted / Nominal";
674  std::string fullTitleRat = ";" + xlabel + ";" + ylabelRat;
675 
676  TH1* hRat = (TH1*)h->Clone();
677  TH1* hp1Rat = (TH1*)hp1->Clone();
678  TH1* hm1Rat = (TH1*)hm1->Clone();
679  TH1* hp2Rat = (TH1*)hp2->Clone();
680  TH1* hm2Rat = (TH1*)hm2->Clone();
681 
682  int nBins = h->GetNbinsX();
683  double nom_evts = h->Integral(1, nBins);
684  double evts_diff = 0.;
685  double percent_diff = 100.;
686  double percent_diff_p = 100.;
687  double percent_diff_m = 100.;
688  for(int i = 0; i <= nBins+1; ++i) {
689  hRat->SetBinContent(i, 1.);
690 
691  if(i != 0 && i != nBins+1) {
692  evts_diff += std::max(std::abs(hp1->GetBinContent(i) - h->GetBinContent(i)),
693  std::abs(hm1->GetBinContent(i) - h->GetBinContent(i)));
694  }
695  }
696  percent_diff *= (evts_diff/nom_evts);
697  percent_diff_p *= ((hp1->Integral(1, nBins) - nom_evts)/nom_evts);
698  percent_diff_m *= ((nom_evts - hm1->Integral(1, nBins))/nom_evts);
699 
700  fprintf(text, "%s %.2s %.2s: +/-%6.4f%%, +%6.4f%%, -%6.4f%%\n",
701  strs.fSystS.c_str(), strs.fComponent.c_str(), strs.fDet.c_str(),
702  percent_diff, percent_diff_p, percent_diff_m);
703 
704  hp1Rat->Divide(h);
705  hm1Rat->Divide(h);
706  hp2Rat->Divide(h);
707  hm2Rat->Divide(h);
708 
709  double maxval = h->GetBinContent(h->GetMaximumBin());
710  maxval = std::max(maxval, hp1->GetBinContent(hp1->GetMaximumBin()));
711  maxval = std::max(maxval, hm1->GetBinContent(hm1->GetMaximumBin()));
712  maxval = std::max(maxval, hp2->GetBinContent(hp2->GetMaximumBin()));
713  maxval = std::max(maxval, hm2->GetBinContent(hm2->GetMaximumBin()));
714 
715  double minval = h->GetBinContent(h->GetMinimumBin());
716  minval = std::min(minval, hp1->GetBinContent(hp1->GetMinimumBin()));
717  minval = std::min(minval, hm1->GetBinContent(hm1->GetMinimumBin()));
718  minval = std::min(minval, hp2->GetBinContent(hp2->GetMinimumBin()));
719  minval = std::min(minval, hm2->GetBinContent(hm2->GetMinimumBin()));
720 
721  if(maxval < minval) { std::swap(maxval, minval); }
722 
723  double maxvalRat = hRat->GetBinContent(hRat->GetMaximumBin());
724  maxvalRat = std::max(maxvalRat, hp1Rat->GetBinContent(hp1Rat->GetMaximumBin()));
725  maxvalRat = std::max(maxvalRat, hm1Rat->GetBinContent(hm1Rat->GetMaximumBin()));
726  maxvalRat = std::max(maxvalRat, hp2Rat->GetBinContent(hp2Rat->GetMaximumBin()));
727  maxvalRat = std::max(maxvalRat, hm2Rat->GetBinContent(hm2Rat->GetMaximumBin()));
728 
729  double minvalRat = hRat->GetBinContent(hRat->GetMinimumBin());
730  for(int i = 1; i <= nBins; ++i) {
731  if(hp1Rat->GetBinContent(i) > 0.) {
732  minvalRat = std::min(minvalRat, hp1Rat->GetBinContent(i));
733  }
734  if(hm1Rat->GetBinContent(i) > 0.) {
735  minvalRat = std::min(minvalRat, hm1Rat->GetBinContent(i));
736  }
737  if(hp2Rat->GetBinContent(i) > 0.) {
738  minvalRat = std::min(minvalRat, hp2Rat->GetBinContent(i));
739  }
740  if(hm2Rat->GetBinContent(i) > 0.) {
741  minvalRat = std::min(minvalRat, hm2Rat->GetBinContent(i));
742  }
743  }
744 
745  if(std::abs(maxvalRat - 1.) > std::abs(1. - minvalRat)) {
746  minvalRat = 1 - std::abs(maxvalRat - 1.);
747  maxvalRat = 1 + std::abs(maxvalRat - 1.);
748  }
749  else {
750  minvalRat = 1 - std::abs(1. - minvalRat);
751  maxvalRat = 1 + std::abs(1. - minvalRat);
752  }
753 
754  CenterTitles(h);
755  h->SetLineColor(kBlack);
756  h->SetMaximum(1.1*maxval);
757  h->SetMinimum(0);
758  h->SetTitle(fullTitle.c_str());
759 
760  CenterTitles(hp1);
761  hp1->SetLineColor(kRed);
762  hp1->SetLineStyle(2);
763  hp1->SetMaximum(1.1*maxval);
764  hp1->SetMinimum(0);
765  hp1->SetTitle(fullTitle.c_str());
766 
767  CenterTitles(hm1);
768  hm1->SetLineColor(kRed);
769  hm1->SetLineStyle(2);
770  hm1->SetMaximum(1.1*maxval);
771  hm1->SetMinimum(0);
772  hm1->SetTitle(fullTitle.c_str());
773 
774  CenterTitles(hp2);
775  hp2->SetLineColor(kBlue);
776  hp2->SetLineStyle(2);
777  hp2->SetMaximum(1.1*maxval);
778  hp2->SetMinimum(0);
779  hp2->SetTitle(fullTitle.c_str());
780 
781  CenterTitles(hm2);
782  hm2->SetLineColor(kBlue);
783  hm2->SetLineStyle(2);
784  hm2->SetMaximum(1.1*maxval);
785  hm2->SetMinimum(0);
786  hm2->SetTitle(fullTitle.c_str());
787 
788  CenterTitles(hRat);
789  hRat->SetLineColor(kBlack);
790  hRat->SetMaximum(1.1*maxvalRat);
791  hRat->SetMinimum(0.9*minvalRat);
792  hRat->SetTitle(fullTitleRat.c_str());
793 
794  CenterTitles(hp1Rat);
795  hp1Rat->SetLineColor(kRed);
796  hp1Rat->SetLineStyle(2);
797  hp1Rat->SetMaximum(1.1*maxvalRat);
798  hp1Rat->SetMinimum(0.9*minvalRat);
799  hp1Rat->SetTitle(fullTitleRat.c_str());
800 
801  CenterTitles(hm1Rat);
802  hm1Rat->SetLineColor(kRed);
803  hm1Rat->SetLineStyle(2);
804  hm1Rat->SetMaximum(1.1*maxvalRat);
805  hm1Rat->SetMinimum(0.9*minvalRat);
806  hm1Rat->SetTitle(fullTitleRat.c_str());
807 
808  CenterTitles(hp2Rat);
809  hp2Rat->SetLineColor(kBlue);
810  hp2Rat->SetLineStyle(2);
811  hp2Rat->SetMaximum(1.1*maxvalRat);
812  hp2Rat->SetMinimum(0.9*minvalRat);
813  hp2Rat->SetTitle(fullTitleRat.c_str());
814 
815  CenterTitles(hm2Rat);
816  hm2Rat->SetLineColor(kBlue);
817  hm2Rat->SetLineStyle(2);
818  hm2Rat->SetMaximum(1.1*maxvalRat);
819  hm2Rat->SetMinimum(0.9*minvalRat);
820  hm2Rat->SetTitle(fullTitleRat.c_str());
821 
822  double xL = 0.6, xR = 0.85;
823  double yB = 0.6, yT = 0.85;
824  TLegend* leg = new TLegend(xL, yB, xR, yT);
825  leg->AddEntry(h, "Nominal", "l");
826  leg->AddEntry(hp1, "+/- 1 #sigma", "l");
827  leg->AddEntry(hp2, "+/- 2 #sigma", "l");
828 
829  TCanvas* c = new TCanvas(name.c_str(), title.c_str(), 800, 800);
830  TPad* pSpecs = new TPad("pSpecs", "", 0., 0.375, 1., 1.);
831  TPad* pRatio = new TPad("pRatio", "", 0., 0., 1., 0.375);
832  pSpecs->Draw();
833  pRatio->Draw();
834 
835  pSpecs->cd();
836  h->Draw("hist");
837  hp2->Draw("hist same");
838  hm2->Draw("hist same");
839  hp1->Draw("hist same");
840  hm1->Draw("hist same");
841  leg->Draw();
842  Simulation();
843 
844  pRatio->cd();
845  hRat->Draw("hist");
846  hp2Rat->Draw("hist same");
847  hm2Rat->Draw("hist same");
848  hp1Rat->Draw("hist same");
849  hm1Rat->Draw("hist same");
850 
851  c->Update();
852  out->WriteTObject(c);
853  return;
854  }
855 
856  //----------------------------------------------------------------------------
857  /// Plot ND syst shifted and nominal spectra
858  /// This is for plotting an on/off style shift
859  void PlotMultiSyst(TH1* nom, std::vector<IDecomp*> shifts,
860  TDirectory* out, FILE* text, strings strs,
861  double POT, bool NC, bool use_max)
862  {
863 
864  std::vector<TH1*> hp;
865  for( auto shift: shifts){
866  TH1* hp1 =(NC ? GetNC(shift, POT) : GetBG(shift, POT));
867  hp.push_back(hp1);
868 }
869  // AddErrorInQuadrature(nom, hp, use_max);
870  PlotMultiSyst(nom, hp, out, text, strs);
871  return;
872  }
873 
874  //----------------------------------------------------------------------------
875  /// Plot FD syst shifted and nominal spectra
876  /// Should be able to change the input to IPrediction
877  /// This is for plotting an on/off style shift
878  void PlotMultiSyst(TH1* nom, std::vector<PredictionNoExtrap*> shifts,
879  TDirectory* out, FILE* text, strings strs,
880  double POT, bool NC, bool use_max)
881  {
882  std::vector<TH1*> hp;
883 
884 for(auto shift: shifts){
885  TH1* hp1=(NC ? GetNC(shift, POT) : GetBG(shift, POT));
886 
887  hp.push_back(hp1);
888 }
889 // AddErrorInQuadrature(nom, hp, use_max);
890  PlotMultiSyst(nom, hp, out, text, strs);
891  return;
892  }
893 
894  //----------------------------------------------------------------------------
895  /// Plot FD syst shifted and nominal spectra
896  /// Should be able to change the input to IPrediction and make this redundant
897  /// This is for plotting an on/off style shift
898  void PlotMultiSyst(TH1* nom, std::vector<PredictionSterile*> shifts,
899  TDirectory* out, FILE* text, strings strs,
900  double POT, bool NC, bool use_max)
901  {
902  std::vector<TH1*> hp;
903  for(auto shift: shifts){
904  TH1* hp1 =(NC ? GetNC(shift, POT) : GetBG(shift, POT));
905  hp.push_back(hp1);
906 }
907  // AddErrorInQuadrature(nom, hp, use_max);
908  PlotMultiSyst(nom, hp, out, text, strs);
909  return;
910  }
911 
912  //----------------------------------------------------------------------------
913  /// Plot ND syst shifted and nominal spectra
914  /// This is for plotting a syst with multiple sigma levels
915  void PlotSyst(TH1* nom, std::map<int, IDecomp*> shifts,
916  TDirectory* out, FILE* text, strings strs,
917  double POT, bool NC, bool use_max)
918  {
919  TH1* hp1 = (NC ? GetNC(shifts[+1], POT) : GetBG(shifts[+1], POT));
920  TH1* hm1 = (NC ? GetNC(shifts[-1], POT) : GetBG(shifts[-1], POT));
921  TH1* hp2 = (NC ? GetNC(shifts[+2], POT) : GetBG(shifts[+2], POT));
922  TH1* hm2 = (NC ? GetNC(shifts[-2], POT) : GetBG(shifts[-2], POT));
923 
924  AddErrorInQuadrature(nom, hp1, hm1, use_max);
925  PlotSyst(nom, hp1, hm1, hp2, hm2, out, text, strs);
926  return;
927  }
928 
929  //----------------------------------------------------------------------------
930  /// Plot FD syst shifted and nominal spectra
931  /// Should be able to change the input to IPrediction
932  /// This is for plotting a syst with multiple sigma levels
933  void PlotSyst(TH1* nom, std::map<int, PredictionNoExtrap*> shifts,
934  TDirectory* out, FILE* text, strings strs,
935  double POT, bool NC, bool use_max)
936  {
937  TH1* hp1 = (NC ? GetNC(shifts[+1], POT) : GetBG(shifts[+1], POT));
938  TH1* hm1 = (NC ? GetNC(shifts[-1], POT) : GetBG(shifts[-1], POT));
939  TH1* hp2 = (NC ? GetNC(shifts[+2], POT) : GetBG(shifts[+2], POT));
940  TH1* hm2 = (NC ? GetNC(shifts[-2], POT) : GetBG(shifts[-2], POT));
941 
942  AddErrorInQuadrature(nom, hp1, hm1, use_max);
943  PlotSyst(nom, hp1, hm1, hp2, hm2, out, text, strs);
944  return;
945  }
946 
947  //----------------------------------------------------------------------------
948  /// Plot FD syst shifted and nominal spectra
949  /// Should be able to change the input to IPrediction and make this redundant
950  /// This is for plotting a syst with multiple sigma levels
951  void PlotSyst(TH1* nom, std::map<int, PredictionSterile*> shifts,
952  TDirectory* out, FILE* text, strings strs,
953  double POT, bool NC, bool use_max)
954  {
955  TH1* hp1 = (NC ? GetNC(shifts[+1], POT) : GetBG(shifts[+1], POT));
956  TH1* hm1 = (NC ? GetNC(shifts[-1], POT) : GetBG(shifts[-1], POT));
957  TH1* hp2 = (NC ? GetNC(shifts[+2], POT) : GetBG(shifts[+2], POT));
958  TH1* hm2 = (NC ? GetNC(shifts[-2], POT) : GetBG(shifts[-2], POT));
959 
960  AddErrorInQuadrature(nom, hp1, hm1, use_max);
961  PlotSyst(nom, hp1, hm1, hp2, hm2, out, text, strs);
962  return;
963  }
964 
965  //----------------------------------------------------------------------------
966  /// Plot the systematic error band around a nominal spectrum
967  /// The errors of the input histogram should be squared
968  void PlotSystBand(TH1* h, TDirectory* out, FILE* text, strings strs)
969  {
970  // Fix the errors by taking the square root!
971  SqrtError(h);
972 
973  std::string xlabel = strs.fXLabel;
974  std::string ylabel = "Events / " + strs.fPOT + " / 0.25 GeV";
975  std::string title = strs.fDet + " " + strs.fComponent +
976  " Error Band for All " + strs.fSystType + " Systematics";
977  std::string fullTitle = ";" + xlabel + ";" + ylabel;
978 // std::string fullTitle = title + ";" + xlabel + ";" + ylabel;
979  std::string name = "c" + strs.fComponent + strs.fDet +
980  strs.fSystType + "Systs" + strs.fSample;
981  std::string ylabelRat = "Shifted / Nominal";
982  std::string fullTitleRat = ";" + xlabel + ";" + ylabelRat;
983 
984  int nBins = h->GetNbinsX();
985  double maxval = 0.;
986  double nom_evts = h->Integral(1, nBins);
987  double evts_diff = 0.;
988  double percent_diff = 100.;
989  for(int i = 1; i <= nBins; ++i) {
990  maxval = std::max(maxval, h->GetBinContent(i) + h->GetBinError(i));
991  evts_diff += std::abs(h->GetBinError(i));
992  }
993  percent_diff *= (evts_diff/nom_evts);
994 
995  fprintf(text, "\nAll systs %.2s %.2s: +/-%6.4f%%\n",
996  strs.fComponent.c_str(), strs.fDet.c_str(),
997  percent_diff);
998 
999  CenterTitles(h);
1000  h->SetLineColor(kTotalMCColor);
1001  h->SetMaximum(1.1*maxval);
1002  h->SetMinimum(0);
1003  h->SetTitle(fullTitle.c_str());
1004 
1005  TH1* hErr = (TH1*)h->Clone();
1006  CenterTitles(hErr);
1007  hErr->SetFillColor(kTotalMCColor-9);
1008  hErr->SetMaximum(1.1*maxval);
1009  hErr->SetMinimum(0);
1010 
1011  double maxvalRat = 1.;
1012  TH1* hRat = (TH1*)h->Clone();
1013  TH1* hUpp = (TH1*)h->Clone();
1014  TH1* hDwn = (TH1*)h->Clone();
1015  for(int i = 1, n = hRat->GetNbinsX(); i <= n; ++i) {
1016  double binVal = hRat->GetBinContent(i);
1017  double errVal = hRat->GetBinError(i);
1018  if(binVal <= 0.) {
1019  binVal = 1.;
1020  errVal = 0.;
1021  }
1022 
1023  hRat->SetBinContent(i, 1);
1024  hUpp->SetBinContent(i, (binVal + errVal)/binVal);
1025  hDwn->SetBinContent(i, (binVal - errVal)/binVal);
1026 
1027  double diff = (binVal + errVal)/binVal;
1028  if(diff > maxvalRat) {
1029  maxvalRat = diff;
1030  }
1031  }
1032  double minvalRat = 1. - (maxvalRat-1.);
1033 
1034  CenterTitles(hRat);
1035  hRat->SetLineColor(kBlack);
1036  hRat->SetMaximum(1.05*maxvalRat);
1037  hRat->SetMinimum(0.95*minvalRat);
1038  hRat->SetTitle(fullTitleRat.c_str());
1039 
1040  CenterTitles(hUpp);
1041  hUpp->SetLineColor(kRed);
1042  hUpp->SetLineStyle(2);
1043  hUpp->SetMaximum(1.05*maxvalRat);
1044  hUpp->SetMinimum(0.95*minvalRat);
1045  hUpp->SetTitle(fullTitleRat.c_str());
1046 
1047  CenterTitles(hDwn);
1048  hDwn->SetLineColor(kRed);
1049  hDwn->SetLineStyle(2);
1050  hDwn->SetMaximum(1.05*maxvalRat);
1051  hDwn->SetMinimum(0.95*minvalRat);
1052  hDwn->SetTitle(fullTitleRat.c_str());
1053 
1054  TCanvas* c = new TCanvas(name.c_str(), title.c_str(), 800, 800);
1055  TPad* pSpecs = new TPad("pSpecs", "", 0., 0.375, 1., 1.);
1056  TPad* pRatio = new TPad("pRatio", "", 0., 0., 1., 0.375);
1057  pSpecs->Draw();
1058  pRatio->Draw();
1059 
1060  pSpecs->cd();
1061  hErr->Draw("e2");
1062  h->Draw("hist same");
1063  Simulation();
1064 
1065  pRatio->cd();
1066  hRat->Draw("hist");
1067  hUpp->Draw("hist same");
1068  hDwn->Draw("hist same");
1069 
1070  c->Update();
1071  out->WriteTObject(c);
1072  return;
1073  }
1074 
1075  //----------------------------------------------------------------------------
1076  /// Save all of the objects in the input maps to the out directory/file
1077  void SaveMaps(
1078  TDirectory* out,
1079  std::map<std::string, IDecomp*> decomps_nominal,
1080  std::map<std::string, std::map<std::string, std::map<int, IDecomp*> > > decomps_shifted,
1081  std::map<std::string, PredictionNoExtrap*> predsNE_nominal,
1082  std::map<std::string, std::map<std::string, std::map<int, PredictionNoExtrap*> > > predsNE_shifted,
1083  std::map<std::string, PredictionSterile*> predsSt_nominal,
1084  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > > predsSt_shifted
1085  )
1086  {
1087  std::map<std::string, TDirectory*> savedirs;
1088  std::string dir;
1089  std::string sep = "__";
1090 
1091  for(const auto& decomp : decomps_nominal) {
1092  dir = "decomp" + sep + decomp.first; // sample name
1093  savedirs[dir] = out->mkdir(dir.c_str());
1094  decomp.second->SaveTo(savedirs[dir]);
1095  }
1096 
1097  for(const auto& sample : decomps_shifted) {
1098  for(const auto& shiftlabel : sample.second) {
1099  for(const auto& decomp : shiftlabel.second) {
1100  dir = "decomp" + sep + sample.first +
1101  sep + shiftlabel.first + sep +
1102  StringFromInt(decomp.first, false);
1103  savedirs[dir] = out->mkdir(dir.c_str());
1104  decomp.second->SaveTo(savedirs[dir]);
1105  }
1106  }
1107  }
1108 
1109  for(const auto& pred : predsNE_nominal) {
1110  dir = "predNE" + sep + pred.first;
1111  savedirs[dir] = out->mkdir(dir.c_str());
1112  pred.second->SaveTo(savedirs[dir]);
1113  }
1114 
1115  for(const auto& sample : predsNE_shifted) {
1116  for(const auto& shiftlabel : sample.second) {
1117  for(const auto& pred : shiftlabel.second) {
1118  dir = "predNE" + sep + sample.first +
1119  sep + shiftlabel.first + sep +
1120  StringFromInt(pred.first, false);
1121  savedirs[dir] = out->mkdir(dir.c_str());
1122  pred.second->SaveTo(savedirs[dir]);
1123  }
1124  }
1125  }
1126 
1127  for(const auto& pred : predsSt_nominal) {
1128  dir = "predSt" + sep + pred.first;
1129  savedirs[dir] = out->mkdir(dir.c_str());
1130  pred.second->SaveTo(savedirs[dir]);
1131  }
1132 
1133  for(const auto& sample : predsSt_shifted) {
1134  for(const auto& shiftlabel : sample.second) {
1135  for(const auto& pred : shiftlabel.second) {
1136  dir = "predSt" + sep + sample.first +
1137  sep + shiftlabel.first + sep +
1138  StringFromInt(pred.first, false);
1139  savedirs[dir] = out->mkdir(dir.c_str());
1140  pred.second->SaveTo(savedirs[dir]);
1141  }
1142  }
1143  }
1144 
1145  return;
1146  }
1147 }
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.
osc::IOscCalculatorAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
void PlotMultiSyst(TH1 *h, std::vector< TH1 * > hp, TDirectory *out, FILE *text, strings strs)
Plot multiple systematics.
Definition: PPFXHelper.h:464
virtual Spectrum AntiNueComponent() const =0
set< int >::iterator it
std::string fComponent
Definition: PPFXHelper.h:103
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:553
virtual Spectrum PredictComponent(osc::IOscCalculator *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
virtual Spectrum NumuComponent() const =0
Optimized version of OscCalculatorPMNS.
Definition: StanTypedefs.h:28
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
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
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< PredictionSterile > LoadFrom(TDirectory *dir)
std::string fPOT
Definition: PPFXHelper.h:105
virtual Spectrum AntiNumuComponent() const =0
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
::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
virtual Spectrum NCTotalComponent() const
Definition: IDecomp.h:18
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
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
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
static std::unique_ptr< PredictionExtrap > LoadFrom(TDirectory *dir)
TH1 * GetBG(IDecomp *specs, double POT)
Definition: PPFXHelper.h:236
float count
Definition: make_cosmics.C:27
std::string StringFromDouble(double pot)
Definition: PPFXHelper.h:128
virtual Spectrum NueComponent() const =0
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
std::unique_ptr< IDecomp > LoadFrom< IDecomp >(TDirectory *dir)
Definition: IDecomp.cxx:32