CovMxSysts.cxx
Go to the documentation of this file.
1 
3 
4 #include "TKey.h"
5 #include "TVectorD.h"
6 
7 #include <iostream>
8 #include <sstream>
9 
10 namespace ana {
11 
12  using std::cout;
13  using std::endl;
14  using std::make_pair;
15  using std::make_unique;
16  using std::map;
17  using std::ostringstream;
18  using std::pair;
19  using std::runtime_error;
20  using std::string;
21  using std::unique_ptr;
22 
23  // -------------------------------------------------------------------------------------
24  NuISyst::NuISyst(string name, string label, string sampleName,
25  map<int, pair<TH1*, TH1*>> shifts) :
26  ISyst(name, label),
27  fSampleName(sampleName),
28  fShifts(shifts)
29  {}
30 
31  // -------------------------------------------------------------------------------------
33 
34  // Clean up the weight histograms
35  for (auto shift : fShifts) {
36  if (shift.second.first) delete shift.second.first;
37  if (shift.second.second) delete shift.second.second;
38  }
39 
40  } // NuISyst destructor
41 
42  // -------------------------------------------------------------------------------------
43  /// Return base name of systematic
44  string NuISyst::BaseName() const {
45  if (fSampleName.empty()) return ShortName(); // If there's no sample prefix to remove
46  return ShortName().substr(fSampleName.size()+1);
47  }
48 
49  // -------------------------------------------------------------------------------------
50  /// SaveTo implementation for NuISyst
51  void NuISyst::SaveTo(TDirectory* dir) const {
52 
53  TDirectory* tmp = gDirectory;
54  dir->cd();
55  TObjString(ShortName().c_str()).Write("name");
56  TObjString(LatexName().c_str()).Write("label");
57  TObjString(fSampleName.c_str()).Write("samplename");
58  TDirectory* sigmasDir = dir->mkdir("Sigmas");
59 
60  // Save each shift
61  for (auto const& sigma : fShifts) {
62 
63  TDirectory* sigmaDir = sigmasDir->mkdir(Form("Sigma%d", sigma.first));
64 
65  // Save sigma value
66  TVectorD vSigma(1);
67  vSigma[0] = sigma.first;
68  sigmaDir->WriteTObject(&vSigma, "sigma");
69 
70  // Save shift histograms
71  sigmaDir->WriteTObject(sigma.second.first, "sig");
72  sigmaDir->WriteTObject(sigma.second.second, "bkg");
73  }
74  tmp->cd();
75 
76  } // function NuISyst::SaveTo
77 
78  // -------------------------------------------------------------------------------------
79  /// Get weight for a given energy, channel and sigma shift
80  double NuISyst::WeightFor(SystChannel channel, double sigma, double calE) const {
81 
82  // If it's a one-sided shift (ie. no -1 sigma weights) then we want to return
83  // 1 for a negative sigma
84  if (sigma < 0 and not fShifts.count(-1)) return 1;
85  const int bin = fShifts.at(0).first->FindBin(calE);
86 
87  auto low = fShifts.begin();
88  auto end = fShifts.end();
89  std::advance(end, -2);
90  if (sigma < fShifts.begin()->first)
91  low = fShifts.begin();
92  else if (sigma >= end->first)
93  low = end;
94  else {
95  for (auto it = fShifts.begin(); it != end; ++it) {
96  if (sigma >= it->first) {
97  low = it;
98  break;
99  }
100  }
101  }
102 
103  auto high = low;
104  ++high;
105 
106  // // Why would we have templates differing by more than 1 sigma?
107  // // fracpart below assumes this
108  assert(high->first - low->first == 1);
109 
110  TH1* hLow = nullptr;
111  TH1* hHigh = nullptr;
112 
113  if (channel == SystChannel::kSig) {
114  hLow = low->second.first;
115  hHigh = high->second.first;
116  }
117  else {
118  hLow = low->second.second;
119  hHigh = high->second.second;
120  }
121 
122  const double fracpart = sigma - low->first;
123  const double ret = (fracpart*hHigh->GetBinContent(bin)) +
124  ((1-fracpart)*hLow->GetBinContent(bin));
125 
126  return std::max(0., ret); // Keep the LL from blowing up
127 
128  } // function NuISyst::WeightFor
129 
130  // -------------------------------------------------------------------------------------
131  KeySyst::KeySyst(string name, string label)
132  : ISyst(name, label), fOneSided(false)
133  {}
134 
135 } // namespace ana
T max(const caf::Proxy< T > &a, T b)
const XML_Char * name
Definition: expat.h:151
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
set< int >::iterator it
NuISyst(std::string name, std::string label, std::string sample_name="", std::map< int, std::pair< TH1 *, TH1 * > > shifts={})
Definition: CovMxSysts.cxx:24
virtual const std::string & ShortName() const final
The name printed out to the screen.
Definition: ISyst.h:27
virtual double WeightFor(SystChannel channel, double sigma, double calE) const
Get weight for a given energy, channel and sigma shift.
Definition: CovMxSysts.cxx:80
Float_t tmp
Definition: plot.C:36
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
const char * label
SystChannel
Definition: CovMxSysts.h:12
std::string fSampleName
Definition: CovMxSysts.h:31
std::string BaseName() const
Return base name of systematic.
Definition: CovMxSysts.cxx:44
float bin[41]
Definition: plottest35.C:14
double sigma(TH1F *hist, double percentile)
OStream cout
Definition: OStream.cxx:6
virtual const std::string & LatexName() const final
The name used on plots (ROOT&#39;s TLatex syntax)
Definition: ISyst.h:30
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
TDirectory * dir
Definition: macro.C:5
std::map< int, std::pair< TH1 *, TH1 * > > fShifts
Definition: CovMxSysts.h:32
KeySyst(std::string name, std::string label)
Definition: CovMxSysts.cxx:131
assert(nhit_max >=nhit_nbins)
virtual void SaveTo(TDirectory *dir) const
SaveTo implementation for NuISyst.
Definition: CovMxSysts.cxx:51
gm Write()