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, string type) const {
52 
53  cout << "NuISyst::SaveTo implementation is a temporary hack that should be replaced once CommonAna migration is complete." << endl;
54 
55  TDirectory* tmp = gDirectory;
56  dir->cd();
57  TObjString(type.c_str()).Write("type");
58  TObjString(ShortName().c_str()).Write("name");
59  TObjString(LatexName().c_str()).Write("label");
60  TObjString(fSampleName.c_str()).Write("samplename");
61  TDirectory* sigmasDir = dir->mkdir("Sigmas");
62 
63  // Save each shift
64  for (auto const& sigma : fShifts) {
65 
66  TDirectory* sigmaDir = sigmasDir->mkdir(Form("Sigma%d", sigma.first));
67 
68  // Save sigma value
69  TVectorD vSigma(1);
70  vSigma[0] = sigma.first;
71  sigmaDir->WriteTObject(&vSigma, "sigma");
72 
73  // Save shift histograms
74  sigmaDir->WriteTObject(sigma.second.first, "sig");
75  sigmaDir->WriteTObject(sigma.second.second, "bkg");
76  }
77  tmp->cd();
78 
79  } // function NuISyst::SaveTo
80 
81  // -------------------------------------------------------------------------------------
82  /// Get weight for a given energy, channel and sigma shift
83  double NuISyst::WeightFor(SystChannel channel, double sigma, double calE) const {
84 
85  // If it's a one-sided shift (ie. no -1 sigma weights) then we want to return
86  // 1 for a negative sigma
87  if (sigma < 0 and not fShifts.count(-1)) return 1;
88  const int bin = fShifts.at(0).first->FindBin(calE);
89 
90  auto low = fShifts.begin();
91  auto end = fShifts.end();
92  std::advance(end, -2);
93  if (sigma < fShifts.begin()->first)
94  low = fShifts.begin();
95  else if (sigma >= end->first)
96  low = end;
97  else {
98  for (auto it = fShifts.begin(); it != end; ++it) {
99  if (sigma >= it->first) {
100  low = it;
101  break;
102  }
103  }
104  }
105 
106  auto high = low;
107  ++high;
108 
109  // // Why would we have templates differing by more than 1 sigma?
110  // // fracpart below assumes this
111  assert(high->first - low->first == 1);
112 
113  TH1* hLow = nullptr;
114  TH1* hHigh = nullptr;
115 
116  if (channel == SystChannel::kSig) {
117  hLow = low->second.first;
118  hHigh = high->second.first;
119  }
120  else {
121  hLow = low->second.second;
122  hHigh = high->second.second;
123  }
124 
125  const double fracpart = sigma - low->first;
126  const double ret = (fracpart*hHigh->GetBinContent(bin)) +
127  ((1-fracpart)*hLow->GetBinContent(bin));
128 
129  return std::max(0., ret); // Keep the LL from blowing up
130 
131  } // function NuISyst::WeightFor
132 
133  // -------------------------------------------------------------------------------------
134  KeySyst::KeySyst(string name, string label)
135  : ISyst(name, label), fOneSided(false)
136  {}
137 
138  // -------------------------------------------------------------------------------------
139  void covmx::NormSyst::Shift(double sigma, caf::SRProxy* sr, double& weight) const {
140  weight *= 1 + (sigma * fScale);
141  if (weight < 0) weight = 0;
142  }
143 
144 } // namespace ana
145 
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
const Var weight
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
virtual void SaveTo(TDirectory *dir, std::string type) const
SaveTo implementation for NuISyst.
Definition: CovMxSysts.cxx:51
virtual double WeightFor(SystChannel channel, double sigma, double calE) const
Get weight for a given energy, channel and sigma shift.
Definition: CovMxSysts.cxx:83
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
caf::StandardRecord * sr
std::string BaseName() const
Return base name of systematic.
Definition: CovMxSysts.cxx:44
float bin[41]
Definition: plottest35.C:14
virtual void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: CovMxSysts.cxx:139
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
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:134
assert(nhit_max >=nhit_nbins)
gm Write()
enum BeamMode string