CutOptimization.cxx
Go to the documentation of this file.
2 
3 #include "TVectorD.h"
4 #include "TCanvas.h"
5 #include "TH1.h"
6 #include "TH2.h"
7 #include "TH3.h"
8 #include "TFile.h"
9 
10 #include <iostream>
11 
12 namespace ana
13 {
14  void
18  std::string name) {
19  TDirectory * tmp = gDirectory;
20  TCanvas * c = new TCanvas();
21  c->SetLeftMargin(0.15);
22  c->SetRightMargin(0.01);
23 
24 
25  c->cd();
26 
27  hist->Draw(draw_option.c_str());
28  hist->SetTitle(title.c_str());
29  c->Print(name.c_str());
30 
31  tmp->cd();
32  }
33 
34  void
35  CutOptimization::PlotDebug(std::vector<TH1 *> hists,
36  std::vector<std::string> labels,
40  double ymax) {
41  TDirectory * tmp = gDirectory;
42  TCanvas * c = new TCanvas();
43  c->SetLeftMargin(0.15);
44  c->SetRightMargin(0.01);
45 
46  c->cd();
47  int colors[] = {kBlue, kViolet, kMagenta,
48  kPink, kGreen, kTeal,
49  kCyan, kAzure, kGray,
50  kGray+1, kGray+3};
51 
52  TLegend * leg = new TLegend(0.4,0.7,0.7,0.9);
53 
54  double max_val = -1e9;
55  for(auto ihist = 0u; ihist < hists.size(); ihist++) {
56  if(hists[ihist] == NULL) continue;
57  max_val = std::max(max_val, hists[ihist]->GetMaximum());
58  }
59 
60  int drawn = 0;
61  for(auto ihist = 0u; ihist < hists.size(); ihist++) {
62  // pointer could be null
63  if(hists[ihist] == NULL) continue;
64 
65  // if valid pointer, style hist
66  hists[ihist]->SetLineColor(colors[ihist]);
67  hists[ihist]->SetTitle(title.c_str());
68 
69  leg->AddEntry(hists[ihist], labels[ihist].c_str(), "l");
70  // use the axes of first valid histogram
71  if(drawn == 0) {
72  // adjust range so all hists will show up
73  if(ymax < 0)
74  hists[ihist]->GetYaxis()->SetRangeUser(0, max_val * 1.2);
75  else
76  hists[ihist]->GetYaxis()->SetRangeUser(0, ymax);
77  hists[ihist]->Draw(draw_option.c_str());
78  }
79  else {
80  hists[ihist]->Draw((draw_option + " same").c_str());
81  }
82  drawn++;
83  }
84 
85  leg->Draw("same");
86  c->Print(name.c_str());
87 
88  tmp->cd();
89  }
90 
91  void
92  CutOptimization::PlotDebug(const std::vector<TH1 *> universes,
93  Spectrum * nominal,
96  {
97  TH1 * h_nominal = nominal->ToTH1(nominal->POT());
98  PlotDebug(universes, h_nominal, title, name);
99  }
100 
101  void
102  CutOptimization::PlotDebug(const std::vector<TH1 *> universes,
103  TH1 * h_nominal,
106  {
107  TDirectory * tmp = gDirectory;
108  TCanvas * c = new TCanvas();
109  c->SetLeftMargin(0.15);
110  c->SetRightMargin(0.01);
111  c->cd();
112 
113  double max_val = -1e9;
114  for(auto ihist = 0u; ihist < universes.size(); ihist++) {
115  max_val = std::max(max_val, universes[ihist]->GetMaximum());
116  }
117  h_nominal->GetYaxis()->SetRangeUser(0, max_val * 1.2);
118  h_nominal->SetTitle(title.c_str());
119  h_nominal->Draw("hist");
120 
121  for(auto ihist = 0u; ihist < universes.size(); ihist++) {
122  universes[ihist]->Draw("hist same");
123  }
124 
125  h_nominal->Draw("hist same");
126 
127 
128 
129  c->Print(name.c_str());
130 
131  tmp->cd();
132  }
133 
135  const Cut * signal_cut,
136  const Cut * selection_cut)
137  : fSignalCut(signal_cut), fSelectionCut(selection_cut), fOptiHistAxis(opti)
138  {}
139 
140  void
142  {
143  fNominalSystDef = nominal;
144 
150  }
151 
152 
153  void
155  {
156  fSystDefs.push_back(syst);
157 
158  // for one sided systs, store in .up
159  fSignal [syst].up = syst->BuildSpectrumUp(*fSignalCut);
160  fBkgd [syst].up = syst->BuildSpectrumUp(!*fSignalCut);
162  fSelectedBkgd [syst].up = syst->BuildSpectrumUp(*fSelectionCut && !*fSignalCut);
163  fSelected [syst].up = syst->BuildSpectrumUp(*fSelectionCut);
164 
165  if(syst->IsTwoSided()) {
166  fSignal [syst].down = syst->BuildSpectrumDown(*fSignalCut);
167  fBkgd [syst].down = syst->BuildSpectrumDown(!*fSignalCut);
169  fSelectedBkgd [syst].down = syst->BuildSpectrumDown(*fSelectionCut && !*fSignalCut);
170  fSelected [syst].down = syst->BuildSpectrumDown(*fSelectionCut);
171  }
172  }
173 
174  void
176  {
177  for(auto iload = 0u; iload < fNominalSystDef->LoadersUp().size(); iload++)
178  fNominalSystDef->Loaders()[iload]->SetSpillCut(spillcut);
179 
180  for(auto isyst = 0u; isyst < fSystDefs.size(); isyst++) {
181  for(auto iload = 0u; iload < fSystDefs[isyst]->LoadersUp().size(); iload++)
182  fSystDefs[isyst]->LoadersUp()[iload]->SetSpillCut(spillcut);
183  if(fSystDefs[isyst]->IsTwoSided()) {
184  for(auto iload = 0u; iload < fSystDefs[isyst]->LoadersDown().size(); iload++)
185  fSystDefs[isyst]->LoadersDown()[iload]->SetSpillCut(spillcut);
186  }
187  }
188  }
189 
190 
191  //////////////////////////////////////////////////////////////
192 
193  void
195  std::vector<Var> mv_weights)
196  {
197  fMVSystDefs.push_back(syst);
198  fMVSignal[syst] = new Multiverse(syst->Loaders(),
199  *fOptiHistAxis,
200  *fSignalCut,
201  *syst->Shifts(),
202  mv_weights,
203  *syst->Weight());
204 
205  fMVBkgd[syst] = new Multiverse(syst->Loaders(),
206  *fOptiHistAxis,
207  !*fSignalCut,
208  *syst->Shifts(),
209  mv_weights,
210  *syst->Weight());
211  fMVSelectedSignal[syst] = new Multiverse(syst->Loaders(),
212  *fOptiHistAxis,
214  *syst->Shifts(),
215  mv_weights,
216  *syst->Weight());
217  fMVSelectedBkgd[syst] = new Multiverse(syst->Loaders(),
218  *fOptiHistAxis,
220  *syst->Shifts(),
221  mv_weights,
222  *syst->Weight());
223  fMVSelected[syst] = new Multiverse(syst->Loaders(),
224  *fOptiHistAxis,
225  *fSelectionCut,
226  *syst->Shifts(),
227  mv_weights,
228  *syst->Weight());
229  }
230 
231  //////////////////////////////////////////////////////////////
232  void
234  std::vector<SystShifts> mv_shifts) {
235  fMVSystDefs.push_back(syst);
236  fMVSignal[syst] = new Multiverse(syst->Loaders(),
237  *fOptiHistAxis,
238  *fSignalCut,
239  mv_shifts,
240  *syst->Weight());
241  fMVBkgd[syst] = new Multiverse(syst->Loaders(),
242  *fOptiHistAxis,
243  !*fSignalCut,
244  mv_shifts,
245  *syst->Weight());
246  fMVSelectedSignal[syst] = new Multiverse(syst->Loaders(),
247  *fOptiHistAxis,
249  mv_shifts,
250  *syst->Weight());
251  fMVSelectedBkgd[syst] = new Multiverse(syst->Loaders(),
252  *fOptiHistAxis,
254  mv_shifts,
255  *syst->Weight());
256  fMVSelected[syst] = new Multiverse(syst->Loaders(),
257  *fOptiHistAxis,
258  *fSelectionCut,
259  mv_shifts,
260  *syst->Weight());
261  }
262 
263  //////////////////////////////////////////////////////////////
264  void
266  {
267  fNominalSystDef->Go();
268  for(auto isyst = 0u; isyst < fSystDefs.size(); isyst++)
269  fSystDefs[isyst]->Go();
270  for(auto isyst = 0u; isyst < fMVSystDefs.size(); isyst++)
271  fMVSystDefs[isyst]->Go();
272  }
273 
274 
275  //////////////////////////////////////////////////////////////
276  void
277  CutOptimization::SaveTo(TDirectory * dir, const std::string& name) const
278  {
279  TDirectory * tmp = gDirectory;
280 
281  dir = dir->mkdir(name.c_str()); // switch to subdir
282  dir->cd();
283 
284  TObjString("CutOptimization").Write("type");
285 
286  TDirectory * nom_dir = dir->mkdir("Nominal");
287  // save nominal
288  fNominalSignal ->SaveTo(nom_dir, "Signal");
289  fNominalBkgd ->SaveTo(nom_dir, "Bkgd");
290  fNominalSelectedSignal->SaveTo(nom_dir, "SelectedSignal");
291  fNominalSelectedBkgd ->SaveTo(nom_dir, "SelectedBkgd");
292  fNominalSelected ->SaveTo(nom_dir, "Selected");
293 
294  dir->cd();
295 
296  TVectorD nsysts(1);
297  nsysts[0] = fSystDefs.size();
298  nsysts.Write("nsysts");
299 
300  // for each syst def, save up down pairs
301  // Iterating over the vector insures systematics will be saved in the same
302  // order everytime. Critical for hadd_cafana
303  for(auto isyst = 0u; isyst < fSystDefs.size(); isyst++) {
304  // save systdefs
305  fSystDefs[isyst]->SaveTo(dir, TString::Format("syst_defs/%d", isyst).Data());
306 
307  TDirectory * syst_dir = dir->mkdir(TString::Format("syst_%s",fSystDefs[isyst]->GetName().c_str()).Data());
308  SaveToUpDownSpectra( fSignal .at(fSystDefs[isyst]), syst_dir->mkdir("fSignal" ));
309  SaveToUpDownSpectra( fBkgd .at(fSystDefs[isyst]), syst_dir->mkdir("fBkgd" ));
310  SaveToUpDownSpectra( fSelectedSignal.at(fSystDefs[isyst]), syst_dir->mkdir("fSelectedSignal"));
311  SaveToUpDownSpectra( fSelectedBkgd .at(fSystDefs[isyst]), syst_dir->mkdir("fSelectedBkgd" ));
312  SaveToUpDownSpectra( fSelected .at(fSystDefs[isyst]), syst_dir->mkdir("fSelected" ));
313  }
314 
315 
316  TVectorD nmultiverse(1);
317  nsysts[0] = fMVSystDefs.size();
318  nsysts.Write("nmultiverse");
319 
320  // for each mv syst def, save multiverse
321  // Iterating over the vector insures multiverses will be saved in the same
322  // order everytime. Critical for hadd_cafana
323  for(auto imv = 0u; imv < fMVSystDefs.size(); imv++) {
324  // save mv syst_defs
325  fMVSystDefs[imv]->SaveTo(dir, TString::Format("mv_defs/%d", imv).Data());
326 
327  TDirectory * mv_dir = dir->mkdir(TString::Format("mv_%s",fMVSystDefs[imv]->GetName().c_str()).Data());
328  fMVSignal .at(fMVSystDefs[imv])->SaveTo(mv_dir, "fMVSignal" );
329  fMVBkgd .at(fMVSystDefs[imv])->SaveTo(mv_dir, "fMVBkgd" );
330  fMVSelectedSignal.at(fMVSystDefs[imv])->SaveTo(mv_dir, "fMVSelectedSignal");
331  fMVSelectedBkgd .at(fMVSystDefs[imv])->SaveTo(mv_dir, "fMVSelectedBkgd" );
332  fMVSelected .at(fMVSystDefs[imv])->SaveTo(mv_dir, "fMVSelected" );
333  }
334 
335  dir->Write();
336  delete dir;
337 
338  tmp->cd();
339 
340  }
341  //////////////////////////////////////////////////////////////
342  std::unique_ptr<CutOptimization>
344  {
345  dir = dir->GetDirectory(name.c_str()); // switch to subdir
346  assert(dir);
347 
348  Spectrum * nominal_signal;
349  Spectrum * nominal_bkgd;
350  Spectrum * nominal_selected_signal;
351  Spectrum * nominal_selected_bkgd;
352  Spectrum * nominal_selected;
353 
354  std::vector<SystematicDef *> syst_defs;
355  std::map<SystematicDef *, UpDownPair<Spectrum> > syst_signal;
356  std::map<SystematicDef *, UpDownPair<Spectrum> > syst_bkgd;
357  std::map<SystematicDef *, UpDownPair<Spectrum> > syst_selected_signal;
358  std::map<SystematicDef *, UpDownPair<Spectrum> > syst_selected_bkgd;
359  std::map<SystematicDef *, UpDownPair<Spectrum> > syst_selected;
360 
361  std::vector<SystematicDef *> mv_defs;
362  std::map<SystematicDef *, Multiverse * > mv_signal;
363  std::map<SystematicDef *, Multiverse * > mv_bkgd;
364  std::map<SystematicDef *, Multiverse * > mv_selected_signal;
365  std::map<SystematicDef *, Multiverse * > mv_selected_bkgd;
366  std::map<SystematicDef *, Multiverse * > mv_selected;
367 
368  TDirectory * nom_dir = dir->GetDirectory("Nominal");
369  nominal_signal = Spectrum::LoadFrom(nom_dir, "Signal" ).release();
370  nominal_bkgd = Spectrum::LoadFrom(nom_dir, "Bkgd" ).release();
371  nominal_selected_signal = Spectrum::LoadFrom(nom_dir, "SelectedSignal").release();
372  nominal_selected_bkgd = Spectrum::LoadFrom(nom_dir, "SelectedBkgd" ).release();
373  nominal_selected = Spectrum::LoadFrom(nom_dir, "Selected" ).release();
374 
375  TVectorD * vsysts = (TVectorD*) dir->Get("nsysts");
376  int nsysts = (*vsysts)[0];
377  for(int isyst = 0; isyst < nsysts; isyst++) {
378  syst_defs.push_back(SystematicDef::LoadFrom(dir, TString::Format("syst_defs/%d", isyst).Data()).release());
379 
380  TDirectory * syst_dir = dir->GetDirectory(TString::Format("syst_%s", syst_defs.back()->GetName().c_str()));
381  syst_signal [syst_defs.back()] = LoadFromUpDownSpectra(syst_dir->GetDirectory("fSignal" ));
382  syst_bkgd [syst_defs.back()] = LoadFromUpDownSpectra(syst_dir->GetDirectory("fBkgd" ));
383  syst_selected_signal[syst_defs.back()] = LoadFromUpDownSpectra(syst_dir->GetDirectory("fSelectedSignal" ));
384  syst_selected_bkgd [syst_defs.back()] = LoadFromUpDownSpectra(syst_dir->GetDirectory("fSelectedBkgd" ));
385  syst_selected [syst_defs.back()] = LoadFromUpDownSpectra(syst_dir->GetDirectory("fSelected" ));
386  }
387 
388  TVectorD * vmultiverse = (TVectorD*) dir->Get("nmultiverse");
389  int nmultiverse = (*vmultiverse)[0];
390  for(int imv = 0; imv < nmultiverse; imv++) {
391  mv_defs.push_back(SystematicDef::LoadFrom(dir, TString::Format("mv_defs/%d", imv).Data()).release());
392 
393  TDirectory * mv_dir = dir->GetDirectory(TString::Format("mv_%s", mv_defs.back()->GetName().c_str()));
394  mv_signal [mv_defs.back()] = Multiverse::LoadFrom(mv_dir, "fMVSignal" ).release();
395  mv_bkgd [mv_defs.back()] = Multiverse::LoadFrom(mv_dir, "fMVBkgd" ).release();
396  mv_selected_signal[mv_defs.back()] = Multiverse::LoadFrom(mv_dir, "fMVSelectedSignal" ).release();
397  mv_selected_bkgd [mv_defs.back()] = Multiverse::LoadFrom(mv_dir, "fMVSelectedBkgd" ).release();
398  mv_selected [mv_defs.back()] = Multiverse::LoadFrom(mv_dir, "fMVSelected" ).release();
399  }
400 
401  delete dir;
402 
403  return std::make_unique<CutOptimization> (nominal_signal,
404  nominal_bkgd,
405  nominal_selected_signal,
406  nominal_selected_bkgd,
407  nominal_selected,
408  syst_defs,
409  syst_signal,
410  syst_bkgd,
411  syst_selected_signal,
412  syst_selected_bkgd,
413  syst_selected,
414  mv_defs,
415  mv_signal,
416  mv_bkgd,
417  mv_selected_signal,
418  mv_selected_bkgd,
419  mv_selected);
420  }
421 
422  //////////////////////////////////////////////////////////////
424  Spectrum * nominal_bkgd,
425  Spectrum * nominal_selected_signal,
426  Spectrum * nominal_selected_bkgd,
427  Spectrum * nominal_selected,
428  std::vector<SystematicDef*> syst_defs,
429  std::map<SystematicDef*, UpDownPair<Spectrum > > syst_signal,
430  std::map<SystematicDef*, UpDownPair<Spectrum > > syst_bkgd,
431  std::map<SystematicDef*, UpDownPair<Spectrum > > syst_selected_signal,
432  std::map<SystematicDef*, UpDownPair<Spectrum > > syst_selected_bkgd,
433  std::map<SystematicDef*, UpDownPair<Spectrum > > syst_selected,
434  std::vector<SystematicDef*> mv_defs,
435  std::map<SystematicDef*, Multiverse * > mv_signal,
436  std::map<SystematicDef*, Multiverse * > mv_bkgd,
437  std::map<SystematicDef*, Multiverse * > mv_selected_signal,
438  std::map<SystematicDef*, Multiverse * > mv_selected_bkgd,
439  std::map<SystematicDef*, Multiverse * > mv_selected)
440  : fNominalSignal(nominal_signal),
441  fNominalBkgd(nominal_bkgd),
442  fNominalSelectedSignal(nominal_selected_signal),
443  fNominalSelectedBkgd(nominal_selected_bkgd),
444  fNominalSelected(nominal_selected),
445  fSystDefs(syst_defs),
446  fSignal(syst_signal),
447  fBkgd(syst_bkgd),
448  fSelectedSignal(syst_selected_signal),
449  fSelectedBkgd(syst_selected_bkgd),
450  fSelected(syst_selected),
451  fMVSystDefs(mv_defs),
452  fMVSignal(mv_signal),
453  fMVBkgd(mv_bkgd),
454  fMVSelectedSignal(mv_selected_signal),
455  fMVSelectedBkgd(mv_selected_bkgd),
456  fMVSelected(mv_selected)
457  {}
458 
459  //////////////////////////////////////////////////////////////
461  {
462  this->fSignalCut = new Cut(*rhs.fSignalCut);
463  this->fSelectionCut = new Cut(*rhs.fSignalCut);
464  this->fOptiHistAxis = new HistAxis(*rhs.fOptiHistAxis);
465 
467  this->fNominalSignal = new Spectrum(*rhs.fNominalSignal );
468  this->fNominalBkgd = new Spectrum(*rhs.fNominalBkgd );
471  this->fNominalSelected = new Spectrum(*rhs.fNominalSelected );
472  for(auto isyst = 0u; isyst < rhs.fSystDefs.size(); isyst++) {
473  this->fSystDefs.push_back(new SystematicDef(*rhs.fSystDefs[isyst]));
474  this->fSignal [this->fSystDefs.back()] = CopyUpDownSpectrum(rhs.fSignal .at(rhs.fSystDefs[isyst]));
475  this->fBkgd [this->fSystDefs.back()] = CopyUpDownSpectrum(rhs.fBkgd .at(rhs.fSystDefs[isyst]));
476  this->fSelectedSignal[this->fSystDefs.back()] = CopyUpDownSpectrum(rhs.fSelectedSignal.at(rhs.fSystDefs[isyst]));
477  this->fSelectedBkgd [this->fSystDefs.back()] = CopyUpDownSpectrum(rhs.fSelectedBkgd .at(rhs.fSystDefs[isyst]));
478  this->fSelected [this->fSystDefs.back()] = CopyUpDownSpectrum(rhs.fSelected .at(rhs.fSystDefs[isyst]));
479  }
480  for(auto imv = 0u; imv < rhs.fMVSystDefs.size(); imv++) {
481  this->fMVSystDefs.push_back(new SystematicDef(*rhs.fMVSystDefs[imv]));
482  this->fMVSignal [this->fMVSystDefs.back()] = new Multiverse(*rhs.fMVSignal .at(rhs.fMVSystDefs[imv]));
483  this->fMVBkgd [this->fMVSystDefs.back()] = new Multiverse(*rhs.fMVBkgd .at(rhs.fMVSystDefs[imv]));
484  this->fMVSelectedSignal[this->fMVSystDefs.back()] = new Multiverse(*rhs.fMVSelectedSignal.at(rhs.fMVSystDefs[imv]));
485  this->fMVSelectedBkgd [this->fMVSystDefs.back()] = new Multiverse(*rhs.fMVSelectedBkgd .at(rhs.fMVSystDefs[imv]));
486  this->fMVSelected [this->fMVSystDefs.back()] = new Multiverse(*rhs.fMVSelected .at(rhs.fMVSystDefs[imv]));
487  }
488  }
489 
490 
491  //////////////////////////////////////////////////////////////
493  {
494  this->fSignalCut = std::move(rhs.fSignalCut);
495  this->fSelectionCut = std::move(rhs.fSignalCut);
496  this->fOptiHistAxis = std::move(rhs.fOptiHistAxis);
497 
498  this->fNominalSystDef = std::move(rhs.fNominalSystDef );
499  this->fNominalSignal = std::move(rhs.fNominalSignal );
500  this->fNominalBkgd = std::move(rhs.fNominalBkgd );
501  this->fNominalSelectedSignal = std::move(rhs.fNominalSelectedSignal);
502  this->fNominalSelectedBkgd = std::move(rhs.fNominalSelectedBkgd );
503  this->fNominalSelected = std::move(rhs.fNominalSelected );
504 
505  for(auto isyst = 0u; isyst < rhs.fSystDefs.size(); isyst++) {
506  this->fSystDefs.push_back(std::move(rhs.fSystDefs[isyst]));
507  this->fSignal [this->fSystDefs.back()] = CopyUpDownSpectrum(std::move(rhs.fSignal .at(rhs.fSystDefs[isyst])));
508  this->fBkgd [this->fSystDefs.back()] = CopyUpDownSpectrum(std::move(rhs.fBkgd .at(rhs.fSystDefs[isyst])));
509  this->fSelectedSignal[this->fSystDefs.back()] = CopyUpDownSpectrum(std::move(rhs.fSelectedSignal.at(rhs.fSystDefs[isyst])));
510  this->fSelectedBkgd [this->fSystDefs.back()] = CopyUpDownSpectrum(std::move(rhs.fSelectedBkgd .at(rhs.fSystDefs[isyst])));
511  this->fSelected [this->fSystDefs.back()] = CopyUpDownSpectrum(std::move(rhs.fSelected .at(rhs.fSystDefs[isyst])));
512  }
513  for(auto imv = 0u; imv < rhs.fMVSystDefs.size(); imv++) {
514  this->fMVSystDefs.push_back(std::move(rhs.fMVSystDefs[imv]));
515  this->fMVSignal [this->fMVSystDefs.back()] = std::move(rhs.fMVSignal .at(rhs.fMVSystDefs[imv]));
516  this->fMVBkgd [this->fMVSystDefs.back()] = std::move(rhs.fMVBkgd .at(rhs.fMVSystDefs[imv]));
517  this->fMVSelectedSignal[this->fMVSystDefs.back()] = std::move(rhs.fMVSelectedSignal.at(rhs.fMVSystDefs[imv]));
518  this->fMVSelectedBkgd [this->fMVSystDefs.back()] = std::move(rhs.fMVSelectedBkgd .at(rhs.fMVSystDefs[imv]));
519  this->fMVSelected [this->fMVSystDefs.back()] = std::move(rhs.fMVSelected .at(rhs.fMVSystDefs[imv]));
520  }
521 
522  }
523 
524  //////////////////////////////////////////////////////////////
527  {
528  if(this == &rhs) return *this;
529 
530  this->fSignalCut = new Cut(*rhs.fSignalCut);
531  this->fSelectionCut = new Cut(*rhs.fSignalCut);
532  this->fOptiHistAxis = new HistAxis(*rhs.fOptiHistAxis);
533 
535  this->fNominalSignal = new Spectrum(*rhs.fNominalSignal );
536  this->fNominalBkgd = new Spectrum(*rhs.fNominalBkgd );
539  this->fNominalSelected = new Spectrum(*rhs.fNominalSelected );
540  for(auto isyst = 0u; isyst < rhs.fSystDefs.size(); isyst++) {
541  this->fSystDefs.push_back(new SystematicDef(*rhs.fSystDefs[isyst]));
542  this->fSignal [this->fSystDefs.back()] = CopyUpDownSpectrum(rhs.fSignal .at(rhs.fSystDefs[isyst]));
543  this->fBkgd [this->fSystDefs.back()] = CopyUpDownSpectrum(rhs.fBkgd .at(rhs.fSystDefs[isyst]));
544  this->fSelectedSignal[this->fSystDefs.back()] = CopyUpDownSpectrum(rhs.fSelectedSignal.at(rhs.fSystDefs[isyst]));
545  this->fSelectedBkgd [this->fSystDefs.back()] = CopyUpDownSpectrum(rhs.fSelectedBkgd .at(rhs.fSystDefs[isyst]));
546  this->fSelected [this->fSystDefs.back()] = CopyUpDownSpectrum(rhs.fSelected .at(rhs.fSystDefs[isyst]));
547  }
548  for(auto imv = 0u; imv < rhs.fMVSystDefs.size(); imv++) {
549  this->fMVSystDefs.push_back(new SystematicDef(*rhs.fMVSystDefs[imv]));
550  this->fMVSignal [this->fMVSystDefs.back()] = new Multiverse(*rhs.fMVSignal .at(rhs.fMVSystDefs[imv]));
551  this->fMVBkgd [this->fMVSystDefs.back()] = new Multiverse(*rhs.fMVBkgd .at(rhs.fMVSystDefs[imv]));
552  this->fMVSelectedSignal[this->fMVSystDefs.back()] = new Multiverse(*rhs.fMVSelectedSignal.at(rhs.fMVSystDefs[imv]));
553  this->fMVSelectedBkgd [this->fMVSystDefs.back()] = new Multiverse(*rhs.fMVSelectedBkgd .at(rhs.fMVSystDefs[imv]));
554  this->fMVSelected [this->fMVSystDefs.back()] = new Multiverse(*rhs.fMVSelected .at(rhs.fMVSystDefs[imv]));
555  }
556  return *this;
557  }
558 
559  //////////////////////////////////////////////////////////////
562  {
563  this->fSignalCut = std::move(rhs.fSignalCut);
564  this->fSelectionCut = std::move(rhs.fSignalCut);
565  this->fOptiHistAxis = std::move(rhs.fOptiHistAxis);
566 
567  this->fNominalSystDef = std::move(rhs.fNominalSystDef );
568  this->fNominalSignal = std::move(rhs.fNominalSignal );
569  this->fNominalBkgd = std::move(rhs.fNominalBkgd );
570  this->fNominalSelectedSignal = std::move(rhs.fNominalSelectedSignal);
571  this->fNominalSelectedBkgd = std::move(rhs.fNominalSelectedBkgd );
572  this->fNominalSelected = std::move(rhs.fNominalSelected );
573 
574  for(auto isyst = 0u; isyst < rhs.fSystDefs.size(); isyst++) {
575  this->fSystDefs.push_back(std::move(rhs.fSystDefs[isyst]));
576  this->fSignal [this->fSystDefs.back()] = CopyUpDownSpectrum(std::move(rhs.fSignal .at(rhs.fSystDefs[isyst])));
577  this->fBkgd [this->fSystDefs.back()] = CopyUpDownSpectrum(std::move(rhs.fBkgd .at(rhs.fSystDefs[isyst])));
578  this->fSelectedSignal[this->fSystDefs.back()] = CopyUpDownSpectrum(std::move(rhs.fSelectedSignal.at(rhs.fSystDefs[isyst])));
579  this->fSelectedBkgd [this->fSystDefs.back()] = CopyUpDownSpectrum(std::move(rhs.fSelectedBkgd .at(rhs.fSystDefs[isyst])));
580  this->fSelected [this->fSystDefs.back()] = CopyUpDownSpectrum(std::move(rhs.fSelected .at(rhs.fSystDefs[isyst])));
581  }
582  for(auto imv = 0u; imv < rhs.fMVSystDefs.size(); imv++) {
583  this->fMVSystDefs.push_back(std::move(rhs.fMVSystDefs[imv]));
584  this->fMVSignal [this->fMVSystDefs.back()] = std::move(rhs.fMVSignal .at(rhs.fMVSystDefs[imv]));
585  this->fMVBkgd [this->fMVSystDefs.back()] = std::move(rhs.fMVBkgd .at(rhs.fMVSystDefs[imv]));
586  this->fMVSelectedSignal[this->fMVSystDefs.back()] = std::move(rhs.fMVSelectedSignal.at(rhs.fMVSystDefs[imv]));
587  this->fMVSelectedBkgd [this->fMVSystDefs.back()] = std::move(rhs.fMVSelectedBkgd .at(rhs.fMVSystDefs[imv]));
588  this->fMVSelected [this->fMVSystDefs.back()] = std::move(rhs.fMVSelected .at(rhs.fMVSystDefs[imv]));
589  }
590  return *this;
591  }
592 
593  //////////////////////////////////////////////////////////////
594  void
596  {
597  auto ndims = tmp->GetDimension();
598  auto NX = tmp->GetNbinsX();
599  auto NY = tmp->GetNbinsY();
600  auto NZ = tmp->GetNbinsZ();
601 
602  // Useful for a 2-sided cut on a single variable like vertex position and fiducial
603  // constraints.
604  // Fills a 2D histogram where
605  // -- on the x axis is the value of the lower bound cut
606  // -- on the y axis is the value of the upper bound cut
607  // bin contents are the integral of the 1D input histogram between the bounds corresponding to
608  // upper bin edges
609  if(method == kIntegratedInsideBounds) {
610  assert(ndims == 1 && "I don't know how to handle this histogram for kIntegratedInsideBounds");
611 
612  // don't leak clone of 1D
613  delete result;
614 
615  // allow for variable bin width
616  auto array = tmp->GetXaxis()->GetXbins()->GetArray();
617  if(array) {
618  result = new TH2D(UniqueName().c_str(), "",
619  NX, tmp->GetXaxis()->GetXbins()->GetArray(),
620  NX, tmp->GetXaxis()->GetXbins()->GetArray());
621  }
622  else {
623  auto minx = tmp->GetXaxis()->GetBinLowEdge(1);
624  auto maxx = tmp->GetXaxis()->GetBinLowEdge(NX+1);
625  result = new TH2D(UniqueName().c_str(), "",
626  NX, minx, maxx,
627  NX, minx, maxx);
628 
629  }
630  result->GetYaxis()->SetTitle(TString::Format("%s %s", tmp->GetXaxis()->GetTitle(), "(Upper)"));
631  result->GetXaxis()->SetTitle(TString::Format("%s %s", tmp->GetXaxis()->GetTitle(), "(Lower)"));
632 
633  for(auto ibinlower = 1; ibinlower <= NX; ibinlower++) {
634  for(auto ibinupper = ibinlower; ibinupper <=NX; ibinupper++) {
635  double error = 0;
636  double integral = tmp->IntegralAndError(ibinlower, ibinupper, error);
637  result->SetBinContent(ibinlower, ibinupper, integral);
638  result->SetBinContent(ibinlower, ibinupper, integral);
639  }
640  }
641  return;
642  }
643  else {
644  switch(ndims) {
645  case 1:
646  {
647  for(auto ibinx = 1; ibinx <= NX; ibinx++) {
648  auto x_a = method == kIntegratedDown? 0 : ibinx;
649  auto x_b = method == kIntegratedDown? ibinx : NX+1 ;
650  double error = 0;
651  double integral = tmp->IntegralAndError(x_a, x_b, error);
652  result->SetBinContent(ibinx, integral);
653  result->SetBinError (ibinx, error );
654 
655  }
656  break;
657  }
658  case 2:
659  {
660  // this way, the casts are only performed twice
661  TH2 * result_2d = static_cast<TH2*> (result);
662  const TH2 * tmp_2d = static_cast<const TH2*>(tmp );
663  for(auto ibinx = 1; ibinx <= NX; ibinx++) {
664  for(auto ibiny = 1; ibiny <= NY; ibiny++) {
665  auto x_a = method == kIntegratedDown? 0 : ibinx;
666  auto x_b = method == kIntegratedDown? ibinx : NX+1 ;
667  auto y_a = method == kIntegratedDown? 0 : ibiny;
668  auto y_b = method == kIntegratedDown? ibiny : NY+1 ;
669 
670  double error = 0;
671  double integral = tmp_2d->IntegralAndError(x_a, x_b,
672  y_a, y_b, error);
673  result_2d->SetBinContent(ibinx, ibiny, integral);
674  result_2d->SetBinError (ibinx, ibiny, error );
675  }
676  }
677  result = static_cast<TH1*>(result_2d);
678  break;
679  }
680  case 3:
681  {
682  // this way, the casts are only performed twice
683  TH3 * result_3d = static_cast<TH3*> (result);
684  const TH3 * tmp_3d = static_cast<const TH3*>(tmp );
685  for(auto ibinx = 1; ibinx <= NX; ibinx++) {
686  for(auto ibiny = 1; ibiny <= NY; ibiny++) {
687  for(auto ibinz = 1; ibinz <= NZ; ibinz++) {
688  auto x_a = method == kIntegratedDown? 0 : ibinx;
689  auto x_b = method == kIntegratedDown? ibinx : NX+1 ;
690  auto y_a = method == kIntegratedDown? 0 : ibiny;
691  auto y_b = method == kIntegratedDown? ibiny : NY+1 ;
692  auto z_a = method == kIntegratedDown? 0 : ibinz;
693  auto z_b = method == kIntegratedDown? ibinz : NZ+1 ;
694 
695  double error = 0;
696  double integral = tmp_3d->IntegralAndError(x_a, x_b,
697  y_a, y_b,
698  z_a, z_b, error);
699  result_3d->SetBinContent(ibinx, ibiny, ibinz, integral);
700  result_3d->SetBinError (ibinx, ibiny, ibinz, error );
701  }
702  }
703  }
704  result = static_cast<TH1*>(result_3d);
705  break;
706  }
707  }
708  }
709 
710  }
711 
712  //////////////////////////////////////////////////////////////
713  // Multiverse should already be scaled and integrated.
714  // This function only returns the absolute shifts that result from the multiverse
716  CutOptimization::ToUpDownHist(Multiverse * mv, const TH1 * h_nominal)
717  {
718  TH1 * h_up = mv->GetPlusOneSigmaShift (h_nominal);
719  TH1 * h_down = mv->GetMinusOneSigmaShift(h_nominal);
720  // transform into 1 sided shift as the most conservative shift
721  TH1 * ret = (TH1*) h_up->Clone();
722  for(auto ibinx = 1; ibinx <= ret->GetNbinsX(); ibinx++) {
723  for(auto ibiny = 1; ibiny <= ret->GetNbinsY(); ibiny++) {
724  for(auto ibinz = 1; ibinz <= ret->GetNbinsZ(); ibinz++) {
725  double up = h_up ->GetBinContent(ibinx, ibiny, ibinz);
726  double down = h_down ->GetBinContent(ibinx, ibiny, ibinz);
727  double nominal = h_nominal ->GetBinContent(ibinx, ibiny, ibinz);
728  double abs_shift = std::max(std::abs(up - nominal), std::abs(down - nominal));
729  ret->SetBinContent(ibinx, ibiny, ibinz, abs_shift + nominal);
730  }
731  }
732  }
733  return UpDownPair<TH1>{ret};
734  }
735 
736 
737  //////////////////////////////////////////////////////////////
738  TH1 *
739  CutOptimization::ToHist(const Spectrum * spec, OptimizeMethod_t method, double POT)
740  {
741  TH1 * tmp;
742  TH1 * ret;
743 
744  // this will happen for one sided systs
745  if(spec == NULL) return NULL;
746 
747  // if no POT provided, use the spectrum's POT
748  if(POT < 0) POT = spec->POT();
749 
750  // Call different function based on dimension
751  auto ndims = spec->NDimensions();
752  switch(ndims) {
753  case 1:
754  tmp = spec->ToTH1(POT);
755  break;
756  case 2:
757  tmp = spec->ToTH2(POT);
758  break;
759  case 3:
760  tmp = spec->ToTH3(POT);
761  break;
762  default:
763  std::cout << "Spectrum appears to be empty. Something went terribly wrong here" << std::endl;
764  exit(1);
765  }
766 
767  if(method == kIntegratedUp || method == kIntegratedDown || method == kIntegratedInsideBounds) {
768  ret = (TH1*) tmp->Clone();
769 
770  // integrated result gets saved back to ret
771  Integrate(ret, tmp, method);
772  }
773  else {
774  // default is kBinByBin, just return the histogram then
775  ret = tmp;
776  }
777 
778  return ret;
779  }
780 
781  // Efficiency is calculate either bin by bin for each numerator and denominator
782  // or integrated to/from bin in numerator divided by the total number of signal events.
783  // If using an integration method, get total number of signal events from the
784  // first or last bin of integrated denominator.
785  //
786  // Calculates binomial errors on the efficiency
787  // https://home.fnal.gov/~paterno/images/effic.pdf
788  TH1 *
790  {
791  // numerator should already be integrated appropriately
792  TH1 * eff = (TH1*) num->Clone(UniqueName().c_str());
793 
794  // if using an integration method, denom should also be integrated appropriately
795  // and if so, scale numerator by total number of signal events, which is found
796  // in either the first or last bin of denominator
797  if(method == kIntegratedUp || method == kIntegratedDown || method == kIntegratedInsideBounds) {
798  auto NX = denom->GetNbinsX();
799  auto NY = denom->GetNbinsY();
800  auto NZ = denom->GetNbinsZ();
801  double integral = 0;
802  if(method == kIntegratedUp || method == kIntegratedDown) {
803  integral = method == kIntegratedUp ? denom->GetBinContent(1) : denom->GetBinContent(NX, NY, NZ);
804  }
805  else if(method == kIntegratedInsideBounds) {
806  // in this case, the total number of signal events will be
807  // where the upper bound is at a max and lower bound is at a minimum
808  // y axis is upper bound,
809  // x axis is lower bound
810  integral = denom->GetBinContent(1, NY, 1);
811  }
812 
813  eff->Scale(1/integral);
814  // compute binomial errors by hand
815  for(auto ibinx = 1; ibinx <= NX; ibinx++) {
816  for(auto ibiny = 1; ibiny <= NY; ibiny++) {
817  for(auto ibinz = 1; ibinz <= NZ; ibinz++) {
818  double efficiency = eff->GetBinContent(ibinx, ibiny, ibinz);
819  double binom_error = 1/integral * std::sqrt( efficiency * (1 - efficiency / integral) );
820  eff->SetBinError(ibinx, ibiny, ibinz, binom_error);
821  }
822  }
823  }
824  }
825  else {
826  // else just divide numerator with denominator
827  eff->Divide(eff, denom, 1, 1, "B"); // divide with binomial error propagation
828  }
829 
830  eff->SetMaximum(1.2*eff->GetMaximum());
831 
832  return eff;
833  }
834 
835  //////////////////////////////////////////////////////////////
838  {
839  TH1 * up = Efficiency(num.up, denom.up, method);
840  TH1 * down = num.down == NULL? NULL : Efficiency(num.down, denom.down, method);
841  return UpDownPair<TH1>{up, down};
842  }
843 
844  Multiverse *
846  Multiverse * denom,
847  OptimizeMethod_t method)
848  {
849  // calculate uncertainty on efficiency of a multiverse as the
850  // +- 1 sigma band of the efficiency multiverse
851  Multiverse * eff = new Multiverse(*num);
852  if(method == kIntegratedUp || method == kIntegratedDown || method == kIntegratedInsideBounds) {
853  eff->Transform(*denom,
854  [method](TH1 * univ_num, TH1 * univ_denom)
855  {
856  return CutOptimization::Efficiency(univ_num, univ_denom, method);
857  });
858  }
859  else {
860  eff->Divide(*denom, true);
861  }
862  return eff;
863  }
864 
865  //////////////////////////////////////////////////////////////
866  void
868  {
869  mv->Transform([method](TH1 * univ)
870  {
871  TH1 * result = (TH1*) univ->Clone();
872  CutOptimization::Integrate(result, univ, method);
873  return result;
874  });
875  }
876 
877  //////////////////////////////////////////////////////////////
878  void
880  TDirectory * save_dir, std::string plot_dump,
881  bool debug)
882  {
883  double mc_pot = fNominalSignal->POT();
884  double s = data_pot / mc_pot;
885  // transform the data into histograms
886  // and integrate if method == kIntegratedUp || method == kIntegradedDown || method == kIntegratedInsideBounds
887  TH1 * hNominalSignal = ToHist(fNominalSignal , method, mc_pot ); // scaled to mc for eff
888  TH1 * hNominalSelectedSignal = ToHist(fNominalSelectedSignal, method, mc_pot ); // scaled to mc for eff
889  TH1 * hNominalBkgd = ToHist(fNominalBkgd , method, mc_pot ); // scaled to for mc statistics
890  TH1 * hNominalSelectedBkgd = ToHist(fNominalSelectedBkgd , method, mc_pot ); // scaled to for mc statistics
891  TH1 * hNominalSelected = ToHist(fNominalSelected , method, data_pot); // scaled to data
892  TH1 * hNominalEfficiency = Efficiency(hNominalSelectedSignal, hNominalSignal, method); // scaled to mc
893 
894 
895  // scale the multiverses and calculate efficiencies
896  std::map<SystematicDef*, Multiverse *> fMVEfficiency;
897  for(auto mv = fMVSystDefs.begin(); mv != fMVSystDefs.end(); mv++) {
898  fMVSignal [*mv]->Scale(mc_pot / mc_pot); // for posterity
899  fMVSelectedSignal[*mv]->Scale(mc_pot / mc_pot);
900  fMVBkgd [*mv]->Scale(mc_pot / mc_pot);
901  fMVSelectedBkgd [*mv]->Scale(mc_pot / mc_pot);
902  fMVSelected [*mv]->Scale(data_pot / mc_pot);
903 
904  if(method == kIntegratedUp || method == kIntegratedDown || method == kIntegratedInsideBounds) {
905  Integrate(fMVSignal [*mv], method);
906  Integrate(fMVSelectedSignal[*mv], method);
907  Integrate(fMVBkgd [*mv], method);
908  Integrate(fMVSelectedBkgd [*mv], method);
909  Integrate(fMVSelected [*mv], method);
910  }
911 
912  fMVEfficiency[*mv] = Efficiency(fMVSelectedSignal[*mv], fMVSignal[*mv], method);
913  }
914 
915 
916  std::vector<UpDownPair<TH1 > > hSystSignal;
917  std::vector<UpDownPair<TH1 > > hSystBkgd;
918  std::vector<UpDownPair<TH1 > > hSystSelectedSignal;
919  std::vector<UpDownPair<TH1 > > hSystSelectedBkgd;
920  std::vector<UpDownPair<TH1 > > hSystSelected;
921  std::vector<UpDownPair<TH1 > > hSystEfficiency;
922 
923  std::vector<SystematicDef*> systs;
924  for(auto isyst = 0u; isyst < fSystDefs.size(); isyst++) {
925  SystematicDef * syst = fSystDefs[isyst];
926  systs.push_back(syst);
927  hSystSignal .push_back({ToHist(fSignal .at(syst).up, method, mc_pot ), ToHist(fSignal .at(syst).down, method, mc_pot )});
928  hSystSelectedSignal.push_back({ToHist(fSelectedSignal .at(syst).up, method, mc_pot ), ToHist(fSelectedSignal .at(syst).down, method, mc_pot )});
929 
930  hSystBkgd .push_back({ToHist(fBkgd .at(syst).up, method, mc_pot ), ToHist(fBkgd .at(syst).down, method, mc_pot )});
931  hSystSelectedBkgd .push_back({ToHist(fSelectedBkgd .at(syst).up, method, mc_pot ), ToHist(fSelectedBkgd .at(syst).down, method, mc_pot )});
932  hSystSelected .push_back({ToHist(fSelected .at(syst).up, method, data_pot), ToHist(fSelected .at(syst).down, method, data_pot)});
933  hSystEfficiency .push_back(Efficiency(hSystSelectedSignal.back(), hSystSignal.back(), method));
934  }
935  for(auto imv = 0u; imv < fMVSystDefs.size(); imv++) {
936  SystematicDef * mv = fMVSystDefs[imv];
937  systs.push_back(mv);
938 
939  hSystSignal .push_back(ToUpDownHist(fMVSignal [mv], hNominalSignal ));
940  hSystSelectedSignal.push_back(ToUpDownHist(fMVSelectedSignal [mv], hNominalSelectedSignal));
941 
942  hSystBkgd .push_back(ToUpDownHist(fMVBkgd [mv], hNominalBkgd ));
943  hSystSelectedBkgd .push_back(ToUpDownHist(fMVSelectedBkgd [mv], hNominalSelectedBkgd ));
944  hSystSelected .push_back(ToUpDownHist(fMVSelected [mv], hNominalSelected ));
945  hSystEfficiency .push_back(ToUpDownHist(fMVEfficiency [mv], hNominalEfficiency ));
946  }
947 
948  TH1 * hFracUncertaintySelected = (TH1*) hNominalSelected ->Clone(UniqueName().c_str());
949  TH1 * hFracUncertaintyBkgdStat = (TH1*) hNominalSelectedBkgd ->Clone(UniqueName().c_str());
950  TH1 * hFracUncertaintyBkgdSyst = (TH1*) hNominalBkgd ->Clone(UniqueName().c_str());
951  TH1 * hFracUncertaintyEfficiency = (TH1*) hNominalSignal ->Clone(UniqueName().c_str());
952  TH1 * hFracUncertaintyXSec = (TH1*) hNominalSignal ->Clone(UniqueName().c_str());
953 
954  hFracUncertaintySelected ->Scale(0);
955  hFracUncertaintyBkgdStat ->Scale(0);
956  hFracUncertaintyBkgdSyst ->Scale(0);
957  hFracUncertaintyEfficiency ->Scale(0);
958  hFracUncertaintyXSec ->Scale(0);
959 
960  std::vector<TH1*> vabs_uncert_eff_syst (systs.size());
961  std::vector<TH1*> vabs_uncert_bkgd_syst (systs.size());
962  std::vector<TH1*> vfrac_uncert_eff_syst (systs.size());
963  std::vector<TH1*> vfrac_uncert_bkgd_syst(systs.size());
964  TH1* sel_minus_bkgd = (TH1*) hNominalSelected->Clone(UniqueName().c_str());
965  sel_minus_bkgd->Scale(0);
966  for(auto isyst = 0u; isyst < systs.size(); isyst++) {
967  vabs_uncert_eff_syst [isyst] = (TH1*) hNominalSelected->Clone(UniqueName().c_str());
968  vabs_uncert_bkgd_syst [isyst] = (TH1*) hNominalSelected->Clone(UniqueName().c_str());
969  vfrac_uncert_eff_syst [isyst] = (TH1*) hNominalSelected->Clone(UniqueName().c_str());
970  vfrac_uncert_bkgd_syst[isyst] = (TH1*) hNominalSelected->Clone(UniqueName().c_str());
971 
972  if(hNominalSelected->GetDimension() == 1) {
973  vfrac_uncert_bkgd_syst[isyst] ->GetYaxis()->SetTitle("#frac{s#upoint#delta N^{syst}_{bkgd}}{N_{sel} - s#upoint N_{bkgd}}");
974  vfrac_uncert_eff_syst [isyst] ->GetYaxis()->SetTitle("#frac{#delta#varepsilon}{#varepsilon}");
975  vabs_uncert_bkgd_syst [isyst] ->GetYaxis()->SetTitle("s#upoint#delta N^{syst}_{bkgd}");
976  vabs_uncert_eff_syst [isyst] ->GetYaxis()->SetTitle("#delta#varepsilon");
977  }
978 
979  else if(hNominalSelected->GetDimension() == 2) {
980  vfrac_uncert_bkgd_syst[isyst] ->GetZaxis()->SetTitle("#frac{s#upoint#delta N^{syst}_{bkgd}}{N_{sel} - s#upoint N_{bkgd}}");
981  vfrac_uncert_eff_syst [isyst] ->GetZaxis()->SetTitle("#frac{#delta#varepsilon}{#varepsilon}");
982  vabs_uncert_bkgd_syst [isyst] ->GetZaxis()->SetTitle("s#upoint#delta N^{syst}_{bkgd}");
983  vabs_uncert_eff_syst [isyst] ->GetZaxis()->SetTitle("#delta#varepsilon");
984  }
985  vabs_uncert_eff_syst [isyst] ->Scale(0);
986  vabs_uncert_bkgd_syst [isyst] ->Scale(0);
987  vfrac_uncert_eff_syst [isyst] ->Scale(0);
988  vfrac_uncert_bkgd_syst[isyst] ->Scale(0);
989  }
990 
991  auto NX = hNominalSignal->GetNbinsX();
992  auto NY = hNominalSignal->GetNbinsY();
993  auto NZ = hNominalSignal->GetNbinsZ();
994  for(auto ibinx = 1; ibinx <= NX; ibinx++) {
995  for(auto ibiny = 1; ibiny <= NY; ibiny++) {
996  for(auto ibinz = 1; ibinz <= NZ; ibinz++) {
997  double abs_uncert_bkgd_syst = 0;
998  double abs_uncert_eff = 0;
999 
1000  double nom_sel = hNominalSelected ->GetBinContent(ibinx, ibiny, ibinz); // scaled to data
1001  double nom_bkgd = hNominalSelectedBkgd ->GetBinContent(ibinx, ibiny, ibinz); // scaled to mc
1002  double nom_eff = hNominalEfficiency ->GetBinContent(ibinx, ibiny, ibinz);
1003  double nom_signal = nom_sel - s * nom_bkgd;
1004 
1005  sel_minus_bkgd->SetBinContent(ibinx, ibiny, ibinz, nom_signal);
1006  for(auto isyst = 0u; isyst < systs.size(); isyst++) {
1007  double syst_bkgd = 0;
1008  double syst_eff = 0;
1009  // if two sided, take largest shift
1010  if(systs[isyst]->IsTwoSided()) {
1011 
1012  double up_bkgd = hSystSelectedBkgd[isyst].up ->GetBinContent(ibinx, ibiny, ibinz);
1013  double down_bkgd = hSystSelectedBkgd[isyst].down->GetBinContent(ibinx, ibiny, ibinz);
1014 
1015  double up_eff = hSystEfficiency [isyst].up ->GetBinContent(ibinx, ibiny, ibinz);
1016  double down_eff = hSystEfficiency [isyst].down->GetBinContent(ibinx, ibiny, ibinz);
1017 
1018  syst_bkgd = std::max(std::abs(up_bkgd - nom_bkgd), std::abs(down_bkgd - nom_bkgd));
1019  syst_eff = std::max(std::abs(up_eff - nom_eff ), std::abs(down_eff - nom_eff ));
1020  }
1021  else {
1022  syst_bkgd = std::abs(hSystSelectedBkgd[isyst].up ->GetBinContent(ibinx, ibiny, ibinz) - nom_bkgd);
1023  syst_eff = std::abs(hSystEfficiency [isyst].up ->GetBinContent(ibinx, ibiny, ibinz) - nom_eff );
1024  }
1025 
1026  // scale for MC statistics
1027  syst_bkgd = s * syst_bkgd;
1028 
1029  abs_uncert_bkgd_syst += std::pow(syst_bkgd, 2);
1030  abs_uncert_eff += std::pow(syst_eff , 2);
1031 
1032  vabs_uncert_eff_syst [isyst]->SetBinContent(ibinx, ibiny, ibinz, syst_eff );
1033  vabs_uncert_bkgd_syst [isyst]->SetBinContent(ibinx, ibiny, ibinz, syst_bkgd);
1034  if(nom_signal)
1035  vfrac_uncert_bkgd_syst[isyst]->SetBinContent(ibinx, ibiny, ibinz, syst_bkgd / nom_signal);
1036  else
1037  vfrac_uncert_bkgd_syst[isyst]->SetBinContent(ibinx, ibiny, ibinz, 100);
1038  if(nom_eff)
1039  vfrac_uncert_eff_syst [isyst]->SetBinContent(ibinx, ibiny, ibinz, syst_eff / nom_eff);
1040  else
1041  vfrac_uncert_eff_syst [isyst]->SetBinContent(ibinx, ibiny, ibinz, 100);
1042  } // end loop over systematics
1043 
1044  // calculate fractional uncertainty on efficiency
1045  double frac_uncert_sel = 100;
1046  double frac_uncert_bkgd_stat = 100;
1047  double frac_uncert_bkgd_syst = 100;
1048  double frac_uncert_eff = 100;
1049  if(nom_signal) {
1050  frac_uncert_sel = std::sqrt(nom_sel) / nom_signal;
1051  frac_uncert_bkgd_stat = s * std::sqrt(nom_bkgd) / nom_signal;
1052  frac_uncert_bkgd_syst = std::sqrt(abs_uncert_bkgd_syst) / nom_signal;
1053  }
1054  if(nom_eff)
1055  frac_uncert_eff = std::sqrt(abs_uncert_eff) / nom_eff;
1056 
1057  // finally calculate FOM
1058  double frac_uncert_sigma = std::sqrt(std::pow(frac_uncert_sel , 2) +
1059  std::pow(frac_uncert_bkgd_stat, 2) +
1060  std::pow(frac_uncert_bkgd_syst, 2) +
1061  std::pow(frac_uncert_eff , 2) );
1062 
1063 
1064  // fill histograms
1065  hFracUncertaintyXSec ->SetBinContent(ibinx, ibiny, ibinz, frac_uncert_sigma );
1066  hFracUncertaintySelected ->SetBinContent(ibinx, ibiny, ibinz, frac_uncert_sel );
1067  hFracUncertaintyBkgdStat ->SetBinContent(ibinx, ibiny, ibinz, frac_uncert_bkgd_stat);
1068  hFracUncertaintyBkgdSyst ->SetBinContent(ibinx, ibiny, ibinz, frac_uncert_bkgd_syst);
1069  hFracUncertaintyEfficiency ->SetBinContent(ibinx, ibiny, ibinz, frac_uncert_eff );
1070  }// end loop over ibinz
1071  }// end loop over ibiny
1072  }// end loop over ibinx
1073 
1074 
1075  if(debug) {
1076  //////////////////////////////////// DEBUG //////////////////////////////////
1077  // save internal histograms to file for later plotting
1078  TFile * fdebug = new TFile((plot_dump + "/debug.root").c_str(), "recreate");
1079 
1080  hNominalSelectedBkgd->Write("hNominalSelectedBkgd");
1081  hNominalSelected->Write("hNominalSelected");
1082  hNominalSelectedSignal->Write("hNominalSelectedSignal");
1083  hNominalSignal->Write("hNominalSignal");
1084  hNominalBkgd->Write("hNominalBkgd");
1085  hNominalEfficiency->Write("hNominalEfficiency");
1086 
1087  hFracUncertaintyXSec ->Write("hFracUncertaintyXSec" );
1088  hFracUncertaintySelected ->Write("hFracUncertaintySelected" );
1089  hFracUncertaintyBkgdStat ->Write("hFracUncertaintyBkgdStat" );
1090  hFracUncertaintyBkgdSyst ->Write("hFracUncertaintyBkgdSyst" );
1091  hFracUncertaintyEfficiency ->Write("hFracUncertaintyEfficiency" );
1092 
1093 
1094  for(auto isyst = 0u; isyst < systs.size(); isyst++ ) {
1095  hSystSelectedBkgd [isyst].up->Write((systs[isyst]->GetName() + "_hSystSelectedBkgd" + "_UP").c_str());
1096  hSystSelected [isyst].up->Write((systs[isyst]->GetName() + "_hSystSelected" + "_UP").c_str());
1097  hSystSelectedSignal[isyst].up->Write((systs[isyst]->GetName() + "_hSystSelectedSignal" + "_UP").c_str());
1098  hSystBkgd [isyst].up->Write((systs[isyst]->GetName() + "_hSystBkgd" + "_UP").c_str());
1099  hSystSignal [isyst].up->Write((systs[isyst]->GetName() + "_hSystSignal" + "_UP").c_str());
1100  hSystEfficiency [isyst].up->Write((systs[isyst]->GetName() + "_hSystEfficiency" + "_UP").c_str());
1101 
1102  if(hSystSelectedBkgd[isyst].down != NULL) {
1103  hSystSelectedBkgd [isyst].down->Write((systs[isyst]->GetName() + "_hSystSelectedBkgd" + "_DOWN").c_str());
1104  hSystSelected [isyst].down->Write((systs[isyst]->GetName() + "_hSystSelected" + "_DOWN").c_str());
1105  hSystSelectedSignal[isyst].down->Write((systs[isyst]->GetName() + "_hSystSelectedSignal" + "_DOWN").c_str());
1106  hSystBkgd [isyst].down->Write((systs[isyst]->GetName() + "_hSystBkgd" + "_DOWN").c_str());
1107  hSystSignal [isyst].down->Write((systs[isyst]->GetName() + "_hSystSignal" + "_DOWN").c_str());
1108  hSystEfficiency [isyst].down->Write((systs[isyst]->GetName() + "_hSystEfficiency" + "_DOWN").c_str());
1109  }
1110 
1111  vabs_uncert_eff_syst [isyst]->Write((systs[isyst]->GetName() + "_abs_uncert_eff_syst" ).c_str());
1112  vabs_uncert_bkgd_syst [isyst]->Write((systs[isyst]->GetName() + "_abs_uncert_bkgd_syst" ).c_str());
1113  vfrac_uncert_eff_syst [isyst]->Write((systs[isyst]->GetName() + "_frac_uncert_eff_syst" ).c_str());
1114  vfrac_uncert_bkgd_syst[isyst]->Write((systs[isyst]->GetName() + "_frac_uncert_bkgd_syst").c_str());
1115 
1116 
1117 // debug_syst_sel_bkgd_up .push_back(hSystSelectedBkgd [isyst].up);
1118 // debug_syst_sel_bkgd_down .push_back(hSystSelectedBkgd [isyst].down);
1119 // debug_syst_sel_up .push_back(hSystSelected [isyst].up);
1120 // debug_syst_sel_down .push_back(hSystSelected [isyst].down);
1121 // debug_syst_sel_signal_up .push_back(hSystSelectedSignal[isyst].up);
1122 // debug_syst_sel_signal_down .push_back(hSystSelectedSignal[isyst].down);
1123 
1124 // debug_syst_bkgd_up .push_back(hSystBkgd[isyst].up);
1125 // debug_syst_bkgd_down .push_back(hSystBkgd[isyst].down);
1126 
1127  // debug_syst_signal_up .push_back(hSystSignal[isyst].up);
1128  // debug_syst_signal_down .push_back(hSystSignal[isyst].down);
1129 
1130  // debug_syst_eff_up .push_back(hSystEfficiency[isyst].up);
1131  // debug_syst_eff_down .push_back(hSystEfficiency[isyst].down);
1132 
1133  // debug_syst_bkgd_labels.push_back(systs[isyst]->GetName());
1134  // debug_syst_bkgd_labels_no_nominal.push_back(systs[isyst]->GetName());
1135  }
1136 
1137 
1138  /*
1139  PlotDebug(debug_syst_bkgd_up , debug_syst_bkgd_labels, "hist", "Background +1#sigma" , plot_dump + "/debug_syst_bkgd_up.pdf" );
1140  PlotDebug(debug_syst_bkgd_down , debug_syst_bkgd_labels, "hist", "Background -1#sigma" , plot_dump + "/debug_syst_bkgd_down.pdf" );
1141  PlotDebug(debug_syst_sel_up , debug_syst_bkgd_labels, "hist", "Selected +1#sigma" , plot_dump + "/debug_syst_sel_up.pdf" );
1142  PlotDebug(debug_syst_sel_down , debug_syst_bkgd_labels, "hist", "Selected -1#sigma" , plot_dump + "/debug_syst_sel_down.pdf" );
1143  PlotDebug(debug_syst_sel_bkgd_up , debug_syst_bkgd_labels, "hist", "Selected Background +1#sigma", plot_dump + "/debug_syst_sel_bkgd_up.pdf" );
1144  PlotDebug(debug_syst_sel_bkgd_down , debug_syst_bkgd_labels, "hist", "Selected Background -1#sigma", plot_dump + "/debug_syst_sel_bkgd_down.pdf" );
1145  PlotDebug(debug_syst_sel_signal_up , debug_syst_bkgd_labels, "hist", "Selected Signal +1#sigma" , plot_dump + "/debug_syst_sel_signal_up.pdf" );
1146  PlotDebug(debug_syst_sel_signal_down, debug_syst_bkgd_labels, "hist", "Selected Signal -1#sigma" , plot_dump + "/debug_syst_sel_signal_down.pdf");
1147  PlotDebug(debug_syst_signal_up , debug_syst_bkgd_labels, "hist", "Signal +1#sigma" , plot_dump + "/debug_syst_signal_up.pdf" );
1148  PlotDebug(debug_syst_signal_down , debug_syst_bkgd_labels, "hist", "Signal -1#sigma" , plot_dump + "/debug_syst_signal_down.pdf" );
1149  PlotDebug(debug_syst_eff_up , debug_syst_bkgd_labels, "hist", "Efficiency +1#sigma" , plot_dump + "/debug_syst_eff_up.pdf" );
1150  PlotDebug(debug_syst_eff_down , debug_syst_bkgd_labels, "hist", "Efficiency -1#sigma" , plot_dump + "/debug_syst_eff_down.pdf" );
1151 
1152  PlotDebug(vabs_uncert_eff_syst , debug_syst_bkgd_labels_no_nominal, "hist", "Abs. Syst. Uncert Efficiency" , plot_dump + "/debug_syst_abs_uncert_eff.pdf" );
1153  PlotDebug(vabs_uncert_bkgd_syst , debug_syst_bkgd_labels_no_nominal, "hist", "Abs. Syst. Uncert Background" , plot_dump + "/debug_syst_abs_uncert_bkgd.pdf" );
1154  PlotDebug(vfrac_uncert_eff_syst , debug_syst_bkgd_labels_no_nominal, "hist", "Frac. Syst. Uncert Efficiency", plot_dump + "/debug_syst_frac_uncert_eff.pdf" , 2);
1155  PlotDebug(vfrac_uncert_bkgd_syst, debug_syst_bkgd_labels_no_nominal, "hist", "Frac. Syst. Uncert Background", plot_dump + "/debug_syst_frac_uncert_bkgd.pdf" , 2);
1156  */
1157 
1158  for(auto syst : fMVSystDefs) {
1159 
1160  fMVSignal [syst]->SaveTo(fdebug, "MV_" + syst->GetName() + "_Signal" );
1161  fMVBkgd [syst]->SaveTo(fdebug, "MV_" + syst->GetName() + "_Bkgd" );
1162  fMVSelected [syst]->SaveTo(fdebug, "MV_" + syst->GetName() + "_Selected" );
1163  fMVSelectedSignal[syst]->SaveTo(fdebug, "MV_" + syst->GetName() + "_SelectedSignal");
1164  fMVSelectedBkgd [syst]->SaveTo(fdebug, "MV_" + syst->GetName() + "_SelectedBkgd" );
1165  fMVEfficiency [syst]->SaveTo(fdebug, "MV_" + syst->GetName() + "_Efficiency" );
1166 
1167 // const std::vector<TH1*> mv_signal_hists = fMVSignal [syst]->GetUniversesHist(kCyan);
1168 // const std::vector<TH1*> mv_bkgd_hists = fMVBkgd [syst]->GetUniversesHist(kCyan);
1169 // const std::vector<TH1*> mv_sel_hists = fMVSelected [syst]->GetUniversesHist(kCyan);
1170 // const std::vector<TH1*> mv_sel_signal_hists = fMVSelectedSignal[syst]->GetUniversesHist(kCyan);
1171 // const std::vector<TH1*> mv_sel_bkgd_hists = fMVSelectedBkgd [syst]->GetUniversesHist(kCyan);
1172 
1173 // std::string name = syst->GetName();
1174 
1175 // PlotDebug(mv_signal_hists , hNominalSignal , name + " Multiverse: Signal" , plot_dump + "/debug_" + name + "_signal.pdf" );
1176 // PlotDebug(mv_bkgd_hists , hNominalBkgd , name + " Multiverse: Bkgd" , plot_dump + "/debug_" + name + "_bkgd.pdf" );
1177 // PlotDebug(mv_sel_signal_hists, hNominalSelectedSignal, name + " Multiverse: SelectedSignal", plot_dump + "/debug_" + name + "_sel_signal.pdf");
1178 // PlotDebug(mv_sel_bkgd_hists , hNominalSelectedBkgd , name + " Multiverse: SelectedBkgd" , plot_dump + "/debug_" + name + "_sel_bkgd.pdf" );
1179 // PlotDebug(mv_sel_hists , hNominalSelected , name + " Multiverse: Selected" , plot_dump + "/debug_" + name + "_sel.pdf" );
1180  }
1181 
1182  /*
1183  PlotDebug(hNominalSignal , "hist", "hNominalSignal" , plot_dump + "/debug_nom_sig.pdf");
1184  PlotDebug(hNominalSelectedSignal , "hist", "hNominalSelectedSignal", plot_dump + "/debug_nom_sig_sel.pdf");
1185  PlotDebug(hNominalBkgd , "hist", "hNominalBkgd" , plot_dump + "/debug_nom_bkgd.pdf");
1186  PlotDebug(hNominalSelectedBkgd , "hist", "hNominalSelectedBkgd" , plot_dump + "/debug_nom_bkgd_sel.pdf");
1187  PlotDebug(hNominalSelected , "hist", "hNominalSelected" , plot_dump + "/debug_nom_sel.pdf");
1188  PlotDebug(hNominalEfficiency , "hist", "hNominalEfficiency" , plot_dump + "/debug_nom_eff.pdf");
1189 
1190  PlotDebug(sel_minus_bkgd , "hist", "Nominal Selected - Background", plot_dump + "/debug_nom_sel_minus_bkgd.pdf");
1191  */
1192  fdebug->Close();
1193  //////////////////////////////////// END DEBUG //////////////////////////////////
1194  }
1195 
1196  // plot histograms
1197  double maxy = 1;
1198  hFracUncertaintyXSec ->SetTitle("Fractional Uncertainty on Cross-section" );
1199  hFracUncertaintySelected ->SetTitle("Fractional Stat. Uncertainty on Selection" );
1200  hFracUncertaintyBkgdStat ->SetTitle("Fractional Stat. Uncertainty on Background");
1201  hFracUncertaintyBkgdSyst ->SetTitle("Fractional Syst. Uncertainty on Background");
1202  hFracUncertaintyEfficiency ->SetTitle("Fractional Syst. Uncertainty on Efficiency");
1203 
1204  if(fNominalSignal->NDimensions() == 1 && method != kIntegratedInsideBounds) {
1205 
1206  hFracUncertaintyXSec ->GetYaxis()->SetRangeUser(0, maxy);
1207  hFracUncertaintySelected ->GetYaxis()->SetRangeUser(0, maxy);
1208  hFracUncertaintyBkgdStat ->GetYaxis()->SetRangeUser(0, maxy);
1209  hFracUncertaintyBkgdSyst ->GetYaxis()->SetRangeUser(0, maxy);
1210  hFracUncertaintyEfficiency ->GetYaxis()->SetRangeUser(0, maxy);
1211 
1212  hFracUncertaintyXSec ->GetYaxis()->SetTitle("#frac{#delta#sigma}{#sigma}" );
1213  hFracUncertaintySelected ->GetYaxis()->SetTitle("#frac{#sqrt{N_{sel}}}{N_{sel}-s#upoint N_{bkg}}" );
1214  hFracUncertaintyBkgdStat ->GetYaxis()->SetTitle("#frac{s#upoint#sqrt{N_{bkg}}}{N_{sel}-s#upoint N_{bkg}}");
1215  hFracUncertaintyBkgdSyst ->GetYaxis()->SetTitle("#frac{#sqrt{s#upoint#deltaN^{syst}_{bkg}}}{N_{sel}-s#upoint N_{bkg}}");
1216  hFracUncertaintyEfficiency ->GetYaxis()->SetTitle("#frac{#delta#varepsilon}{#varepsilon}" );
1217 
1218  hFracUncertaintyXSec ->GetYaxis()->CenterTitle();
1219  hFracUncertaintySelected ->GetYaxis()->CenterTitle();
1220  hFracUncertaintyBkgdStat ->GetYaxis()->CenterTitle();
1221  hFracUncertaintyBkgdSyst ->GetYaxis()->CenterTitle();
1222  hFracUncertaintyEfficiency ->GetYaxis()->CenterTitle();
1223 
1224  hFracUncertaintyXSec ->SetLineColor(kRed);
1225  hFracUncertaintySelected ->SetLineColor(kRed);
1226  hFracUncertaintyBkgdStat ->SetLineColor(kRed);
1227  hFracUncertaintyBkgdSyst ->SetLineColor(kRed);
1228  hFracUncertaintyEfficiency ->SetLineColor(kRed);
1229 
1230  // plot
1231  TCanvas * canvas = new TCanvas();
1232  canvas->SetLeftMargin(0.15);
1233  canvas->SetRightMargin(0.01);
1234  canvas->cd();
1236  hFracUncertaintyXSec ->Draw("hist");
1237  canvas ->Print(TString::Format("%s/optimize_%s_frac_uncert_xsec.pdf",
1238  plot_dump.c_str(),
1239  label.c_str()));
1240  hFracUncertaintySelected ->Draw("hist");
1241  canvas ->Print(TString::Format("%s/optimize_%s_frac_uncert_sel.pdf",
1242  plot_dump.c_str(),
1243  label.c_str()));
1244  hFracUncertaintyBkgdStat ->Draw("hist");
1245  canvas ->Print(TString::Format("%s/optimize_%s_frac_uncert_bkgd_stat.pdf",
1246  plot_dump.c_str(),
1247  label.c_str()));
1248  hFracUncertaintyBkgdSyst ->Draw("hist");
1249  canvas ->Print(TString::Format("%s/optimize_%s_frac_uncert_bkgd_syst.pdf",
1250  plot_dump.c_str(),
1251  label.c_str()));
1252  hFracUncertaintyEfficiency ->Draw("hist");
1253  canvas ->Print(TString::Format("%s/optimize_%s_frac_uncert_eff.pdf",
1254  plot_dump.c_str(),
1255  label.c_str()));
1256 
1257 
1258  }
1259  else if(fNominalSignal->NDimensions() == 2 || method == kIntegratedInsideBounds) {
1260  hFracUncertaintyXSec ->GetZaxis()->SetRangeUser(0, maxy);
1261  hFracUncertaintySelected ->GetZaxis()->SetRangeUser(0, maxy);
1262  hFracUncertaintyBkgdStat ->GetZaxis()->SetRangeUser(0, maxy);
1263  hFracUncertaintyBkgdSyst ->GetZaxis()->SetRangeUser(0, maxy);
1264  hFracUncertaintyEfficiency ->GetZaxis()->SetRangeUser(0, maxy);
1265 
1266  hFracUncertaintyXSec ->GetZaxis()->SetTitle("#frac{#delta#sigma}{#sigma}" );
1267  hFracUncertaintySelected ->GetZaxis()->SetTitle("#frac{#sqrt{N_{sel}}}{N_{sel}-N_{bkg}}" );
1268  hFracUncertaintyBkgdStat ->GetZaxis()->SetTitle("#frac{#sqrt{N_{bkg}}}{N_{sel}-N_{bkg}}" );
1269  hFracUncertaintyBkgdSyst ->GetZaxis()->SetTitle("#frac{#deltaN^{syst}_{bkg}}{N_{sel}-N_{bkg}}");
1270  hFracUncertaintyEfficiency ->GetZaxis()->SetTitle("#frac{#delta#varepsilon}{#varepsilon}" );
1271 
1272  hFracUncertaintyXSec ->GetYaxis()->CenterTitle();
1273  hFracUncertaintySelected ->GetYaxis()->CenterTitle();
1274  hFracUncertaintyBkgdStat ->GetYaxis()->CenterTitle();
1275  hFracUncertaintyBkgdSyst ->GetYaxis()->CenterTitle();
1276  hFracUncertaintyEfficiency ->GetYaxis()->CenterTitle();
1277 
1278  // plot
1279  TCanvas * canvas = new TCanvas();
1280  canvas->cd();
1281  std::string labelx = fNominalSignal->GetLabels()[0];
1282  std::string labely = method == kIntegratedInsideBounds? labelx : fNominalSignal->GetLabels()[1];
1283  hFracUncertaintyXSec ->Draw("colz");
1284  canvas ->Print(TString::Format("%s/optimize_%s_%s_frac_uncert_xsec.pdf",
1285  plot_dump.c_str(),
1286  labelx.c_str(),
1287  labely.c_str()));
1288  hFracUncertaintySelected ->Draw("colz");
1289  canvas ->Print(TString::Format("%s/optimize_%s_%s_frac_uncert_sel.pdf",
1290  plot_dump.c_str(),
1291  labelx.c_str(),
1292  labely.c_str()));
1293  hFracUncertaintyBkgdStat ->Draw("colz");
1294  canvas ->Print(TString::Format("%s/optimize_%s_%s_frac_uncert_bkgd_stat.pdf",
1295  plot_dump.c_str(),
1296  labelx.c_str(),
1297  labely.c_str()));
1298  hFracUncertaintyBkgdSyst ->Draw("colz");
1299  canvas ->Print(TString::Format("%s/optimize_%s_%s_frac_uncert_bkgd_syst.pdf",
1300  plot_dump.c_str(),
1301  labelx.c_str(),
1302  labely.c_str()));
1303  hFracUncertaintyEfficiency ->Draw("colz");
1304  canvas ->Print(TString::Format("%s/optimize_%s_%s_frac_uncert_eff.pdf",
1305  plot_dump.c_str(),
1306  labelx.c_str(),
1307  labely.c_str()));
1308  }
1309  else {
1310  std::cout << "Don't know how to draw a " << fNominalSignal->NDimensions() << " dimensional histogram.";
1311  exit(1);
1312  }
1313 
1314  // save
1315  save_dir->cd();
1316  hFracUncertaintyXSec ->Write("frac_uncert_xsec" );
1317  hFracUncertaintySelected ->Write("frac_uncert_sel" );
1318  hFracUncertaintyBkgdStat ->Write("frac_uncert_bkgd_stat");
1319  hFracUncertaintyBkgdSyst ->Write("frac_uncert_bkgd_syst");
1320  hFracUncertaintyEfficiency ->Write("frac_uncert_eff" );
1321  }
1322 
1323  //////////////////////////////////////////////////////////////
1324  TH1 *
1325  CutOptimization::AbsUncertainty(const UpDownPair<TH1> syst, const TH1 * nominal) const
1326  {
1327  TH1 * ret = (TH1*) syst.up->Clone();
1328  int NX = syst.up->GetNbinsX();
1329  int NY = syst.up->GetNbinsY();
1330  int NZ = syst.up->GetNbinsZ();
1331  for(auto ibinx = 1; ibinx <= NX; ibinx++) {
1332  for(auto ibiny = 1; ibiny <= NY; ibiny++) {
1333  for(auto ibinz = 1; ibinz <= NZ; ibinz++) {
1334  double nom = nominal->GetBinContent(ibinx, ibiny, ibinz);
1335 
1336  double up = syst.up->GetBinContent(ibinx, ibiny, ibinz);
1337  double down = nom;
1338  if(syst.down) down = syst.down->GetBinContent(ibinx, ibiny, ibinz);
1339 
1340  double error = std::max(std::abs(up - nom), std::abs(down - nom));
1341  ret->SetBinContent(ibinx, ibiny, ibinz, error);
1342  }
1343  }
1344  }
1345 
1346  return ret;
1347  }
1348 
1349  //////////////////////////////////////////////////////////////
1350  TH1 *
1351  CutOptimization::AbsUncertaintySquared(const UpDownPair<TH1> syst, const TH1 * nominal) const
1352  {
1353  TH1 * ret = AbsUncertainty(syst, nominal);
1354  ret->Multiply(ret);
1355  return ret;
1356  }
1357 
1358 
1359  //////////////////////////////////////////////////////////////
1360  void
1361  CutOptimization::Optimize(OptimizeMethod_t method, FOM_t fom, double data_pot,
1362  TDirectory * save_dir, std::string plot_dump,
1363  bool debug)
1364  {
1365  switch(fom) {
1366  case kdSigmaOverSigma: OptimizedSigmaOverSigma(method, data_pot, save_dir, plot_dump, debug);
1367  break;
1368  default:
1369  std::cout << "Your provided FOM is not yet implemented." << std::endl;
1370  exit(0);
1371  }
1372  }
1373 
1374 
1375 
1376 }
T max(const caf::Proxy< T > &a, T b)
const XML_Char * name
Definition: expat.h:151
UpDownPair< TH1 > ToUpDownHist(Multiverse *mv, const TH1 *h_nominal)
TH2 * ToTH2(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 2D to obtain TH2.
Definition: Spectrum.cxx:165
_HistAxis< Var > HistAxis
Definition: HistAxis.h:103
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const SystShifts * Shifts() const
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
Spectrum * GetPlusOneSigmaShift(const Spectrum *)
Definition: Multiverse.cxx:443
static std::unique_ptr< GenericSystematicDef< SRType > > LoadFrom(TDirectory *dir, const std::string &name)
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:148
std::map< SystematicDef *, UpDownPair< Spectrum > > fSignal
std::map< SystematicDef *, UpDownPair< Spectrum > > fSelectedBkgd
double maxy
std::vector< SpectrumLoaderBase * > Loaders() const
Spectrum * BuildSpectrumDown(const _Cut< SRType > cut)
const _Var< SRType > * Weight() const
CutOptimization(const HistAxis *hist_axis, const Cut *signal_cut, const Cut *selection_cut)
const HistAxis * fOptiHistAxis
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
static std::unique_ptr< CutOptimization > LoadFrom(TDirectory *dir, const std::string &name)
const double mc_pot
std::map< SystematicDef *, Multiverse * > fMVBkgd
static TH1 * ToHist(const Spectrum *spec, OptimizeMethod_t method, double POT=-1)
std::map< SystematicDef *, UpDownPair< Spectrum > > fBkgd
Float_t tmp
Definition: plot.C:36
TString hists[nhists]
Definition: bdt_com.C:3
GenericSystematicDef< caf::SRProxy > SystematicDef
Spectrum * GetMinusOneSigmaShift(const Spectrum *)
Definition: Multiverse.cxx:450
float abs(float number)
Definition: d0nt_math.hpp:39
void DefineNominal(SystematicDef *nominal)
UpDownPair< Spectrum > CopyUpDownSpectrum(const UpDownPair< Spectrum > &copy)
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
CutOptimization & operator=(const CutOptimization &rhs)
const char * label
static void PlotDebug(TH1 *hist, std::string draw_option, std::string title, std::string name)
Double_t ymax
Definition: plot.C:25
const XML_Char * s
Definition: expat.h:262
std::string GetName(int i)
Spectrum * BuildSpectrumUp(const _Cut< SRType > cut)
int colors[6]
Definition: tools.h:1
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:546
std::map< SystematicDef *, Multiverse * > fMVSelectedBkgd
void SaveTo(TDirectory *dir, const std::string &name) const
void DefineMultiverseSystematic(SystematicDef *syst, std::vector< Var > mv_weights)
static TH1 * Efficiency(TH1 *num, TH1 *denom, OptimizeMethod_t method)
Spectrum * fNominalSelectedSignal
_Cut< caf::SRProxy > Cut
Representation of a cut (selection) to be applied to a caf::StandardRecord object.
Definition: Cut.h:96
void SetSpillCuts(const SpillCut &)
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:517
unsigned int NDimensions() const
Definition: Spectrum.h:254
static void Integrate(TH1 *&result, const TH1 *tmp, OptimizeMethod_t method)
std::vector< float > Spectrum
Definition: Constants.h:573
std::map< SystematicDef *, UpDownPair< Spectrum > > fSelectedSignal
double POT() const
Definition: Spectrum.h:219
static std::unique_ptr< Multiverse > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Multiverse.cxx:573
OStream cout
Definition: OStream.cxx:6
std::vector< double > POT
void Optimize(OptimizeMethod_t method, FOM_t fom, double data_pot, TDirectory *save_dir, std::string plot_dump=".", bool debug=false)
void OptimizedSigmaOverSigma(OptimizeMethod_t method, double data_pot, TDirectory *save_dir, std::string plot_dump=".", bool debug=false)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
TH3 * ToTH3(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 3D to obtain TH3.
Definition: Spectrum.cxx:218
TDirectory * dir
Definition: macro.C:5
void SaveToUpDownSpectra(const UpDownPair< Spectrum > &save, TDirectory *dir)
void DefineSystematic(SystematicDef *syst)
int num
Definition: f2_nu.C:119
exit(0)
std::map< SystematicDef *, Multiverse * > fMVSelectedSignal
std::map< SystematicDef *, UpDownPair< Spectrum > > fSelected
assert(nhit_max >=nhit_nbins)
std::vector< SystematicDef * > fSystDefs
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
UpDownPair< Spectrum > LoadFromUpDownSpectra(TDirectory *dir)
TH1 * AbsUncertaintySquared(const UpDownPair< TH1 > syst, const TH1 *nominal) const
Spectrum * BuildSpectrum(const _Cut< SRType > cut)
const std::vector< std::string > & GetLabels() const
Definition: Spectrum.h:255
std::map< SystematicDef *, Multiverse * > fMVSelected
void Transform(const Multiverse &rhs, const std::function< ProductFunc_t > &transform)
Definition: Multiverse.cxx:72
void efficiency()
Definition: efficiency.C:58
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
std::vector< SystematicDef * > fMVSystDefs
std::vector< SpectrumLoaderBase * > LoadersUp() const
SystematicDef * fNominalSystDef
void Divide(Multiverse &, bool=false)
Definition: Multiverse.cxx:139
Spectrum * fNominalSelectedBkgd
std::map< SystematicDef *, Multiverse * > fMVSignal
TH1 * AbsUncertainty(const UpDownPair< TH1 > syst, const TH1 *nominal) const