SystShifts.cxx
Go to the documentation of this file.
2 #include "CAFAna/Core/ISyst.h"
3 #include "CAFAna/Core/Registry.h"
4 
5 #include "Utilities/func/MathUtil.h"
7 
8 #include <cassert>
9 #include <memory>
10 #include <iostream>
11 
12 #include "TDirectory.h"
13 #include "TH1.h"
14 #include "TObjString.h"
15 #include "TString.h"
16 #include "TRandom3.h"
17 
18 namespace ana
19 {
20  // Reserve 0 for unshifted
21  int SystShifts::fgNextID = 1;
22 
23  //----------------------------------------------------------------------
25  : fID(fgNextID++)
26  {}
27 
28  //----------------------------------------------------------------------
29  SystShifts::SystShifts(const ISyst* syst, double shift)
30  : fID(fgNextID++)
31  {
32  fSystsDbl.emplace(syst, shift);
33  }
34 
35  //----------------------------------------------------------------------
37  : fID(fgNextID++)
38  {
39  fSystsStan.emplace(syst, shift);
40  // we're always going to maintain a "double" copy
41  // so that when the Stan cache gets invalidated we can still return something usable.
42  fSystsDbl.emplace(syst, util::GetValAs<double>(shift));
43  }
44 
45  //----------------------------------------------------------------------
46  SystShifts::SystShifts(const std::map<const ISyst*, double>& shifts)
47  : fID(fgNextID++)
48  {
49  for(auto it: shifts) fSystsDbl.emplace(it.first, it.second);
50  }
51 
52  //----------------------------------------------------------------------
53  SystShifts::SystShifts(const std::map<const ISyst*, stan::math::var>& shifts)
54  : fID(fgNextID++)
55  {
56  for(auto it: shifts)
57  {
58  fSystsStan.emplace(it.first, it.second);
59  fSystsDbl.emplace(it.first, util::GetValAs<double>(it.second));
60  }
61  }
62 
63  //----------------------------------------------------------------------
64  std::unique_ptr<SystShifts> SystShifts::Copy() const
65  {
66  return std::make_unique<SystShifts>(*this);
67  }
68 
69  //----------------------------------------------------------------------
70  void SystShifts::SetShift(const ISyst* syst, double shift, bool force)
71  {
72  fID = fgNextID++;
73 
74  // if this slot already exists in the Stan systs, and we're not setting the same value,
75  // some shenanigans are going on that we need to figure out. abort.
76  auto itStan = fSystsStan.find(syst);
77  if(itStan != fSystsStan.end() && itStan->second != shift)
78  {
79  std::cerr << "Error Syst '" << syst->ShortName() << " already has a Stan pull set (" << itStan->second << ") "
80  << "and you're trying to set a different double one (" << shift << ")." << std::endl;
81  std::cerr << "You almost certainly didn't mean to do that." << std::endl;
82  std::cerr << "Abort." << std::endl;
83  abort();
84  }
85 
86  fSystsDbl.erase(syst);
87  if(force || shift != 0.) fSystsDbl.emplace(syst, shift);
88  }
89 
90  //----------------------------------------------------------------------
91  void SystShifts::SetShift(const ISyst* syst, stan::math::var shift)
92  {
93  fSystsStan.erase(syst);
94  // note: _always_ put the syst in, even if the value is 0.
95  // autodiff relies on the calculation happening so as to get the gradient
96  fSystsStan.emplace(syst, shift);
97  SetShift(syst, util::GetValAs<double>(shift), true);
98  }
99 
100  //----------------------------------------------------------------------
101  template <>
102  double SystShifts::GetShift(const ISyst* syst) const
103  {
104  assert(syst);
105 
106  auto it = fSystsDbl.find(syst);
107  return (it == fSystsDbl.end()) ? 0 : it->second;
108  }
109 
110  //----------------------------------------------------------------------
111  template <>
113  {
114  assert(syst);
115 
116  auto it = fSystsStan.find(syst);
117  //assert ( (it == fSystsStan.end()) == (fSystsDbl.find(syst) == fSystsDbl.end()) ); // if you're asking for a Stan syst, and it's not there but a double one is, something went wrong
118  if (it == fSystsStan.end())
119  {
120  auto itDbl = fSystsDbl.find(syst);
121  if (itDbl != fSystsDbl.end())
122  {
123  std::cout << "Warning: creating stan::math::var out of double value for syst '" << syst->ShortName() << "'." << std::endl;
124  std::cout << " If you see this repeatedly, it's likely a problem. (A few times during startup is probably harmless.)" << std::endl;
125  fSystsStan[syst] = stan::math::var(itDbl->second);
126  it = fSystsStan.find(syst);
127  }
128  }
129  return (it == fSystsStan.end()) ? 0 : it->second;
130  }
131 
132  //----------------------------------------------------------------------
134  {
135  fID = 0;
136 
137  fSystsDbl.clear();
138  fSystsStan.clear();
139  }
140 
141  //----------------------------------------------------------------------
142  double SystShifts::Penalty() const
143  {
144  double ret = 0;
145  // Systematics are all expressed in terms of sigmas
146  for(auto it: fSystsDbl) ret += it.second * it.second;
147  return ret;
148  }
149 
150  //----------------------------------------------------------------------
152 
153  //----------------------------------------------------------------------
155  double& weight) const
156  {
157  // always fSystsDbl here because this is only used in the event loop, not in fitting
158  // (so the autodiff'd version is not needed)
159  for(auto it: fSystsDbl) it.first->Shift(it.second, sr, weight);
160  }
161 
162  //----------------------------------------------------------------------
164  {
165  // always fSystsDbl here because this is only used in the event loop, not in fitting
166  // (so the autodiff'd version is not needed)
167  for(auto it: fSystsDbl) it.first->TruthShift(it.second, nu, weight);
168  }
169 
170  //----------------------------------------------------------------------
172  {
173  if(IsNominal()) return "nominal";
174 
176  for(auto it: fSystsDbl){
177  if(!ret.empty()) ret += ",";
178  ret += it.first->ShortName() + TString::Format("=%+g", it.second).Data();
179  }
180 
181  return ret;
182  }
183 
184  //----------------------------------------------------------------------
186  {
187  if(IsNominal()) return "Nominal";
188 
190  for(auto it: fSystsDbl){
191  if(!ret.empty()) ret += ", ";
192  ret += it.first->LatexName() + TString::Format(" = %+g", it.second).Data();
193  }
194 
195  return ret;
196  }
197 
198  //----------------------------------------------------------------------
199  std::vector<const ISyst*> SystShifts::ActiveSysts() const
200  {
201  std::vector<const ISyst*> ret;
202  for(auto it: fSystsDbl) ret.push_back(it.first);
203  return ret;
204  }
205 
206  //----------------------------------------------------------------------
207  void SystShifts::SaveTo(TDirectory* dir) const
208  {
209  TDirectory* tmp = gDirectory;
210 
211  dir->cd();
212  TObjString("SystShifts").Write("type");
213 
214  // Don't write any histogram for the nominal case
215  if(!fSystsDbl.empty()){
216  TH1D h("", "", fSystsDbl.size(), 0, fSystsDbl.size());
217  int ibin = 0;
218  for(auto it: fSystsDbl){
219  ++ibin;
220  h.GetXaxis()->SetBinLabel(ibin, it.first->ShortName().c_str());
221  h.SetBinContent(ibin, it.second);
222  }
223  h.Write("vals");
224  }
225 
226  tmp->cd();
227  }
228 
229  //----------------------------------------------------------------------
230  std::unique_ptr<SystShifts> SystShifts::LoadFrom(TDirectory* dir)
231  {
232  TObjString* tag = (TObjString*)dir->Get("type");
233  assert(tag);
234  assert(tag->GetString() == "SystShifts");
235 
236  auto ret = std::make_unique<SystShifts>();
237 
238  TH1* h = (TH1*)dir->Get("vals");
239  if(h){ // no histogram means nominal
240  for(int i = 1; i <= h->GetNbinsX(); ++i){
241  ret->SetShift(Registry<ISyst>::ShortNameToPtr(h->GetXaxis()->GetBinLabel(i)),
242  h->GetBinContent(i));
243  }
244  }
245 
246  return ret;
247  }
248 
249  //----------------------------------------------------------------------
250  std::vector <SystShifts> GetSystShiftsMultiverse(
251  const std::vector <std::pair <const ISyst*, SystMode> > systConfigs,
252  const int nUniverses, const int seed)
253  {
254  std::vector <SystShifts> shifts = {};
255  int tempseed = seed;
256 
257  for (int univ=0; univ<nUniverses; ++univ){
258  SystShifts tempshift;
259 
260  for (auto systConfig:systConfigs) {
261  TRandom3 rand (tempseed++);
262  if(systConfig.second == kSystExclude) continue;
263  else if (systConfig.second == kSystGaussian) {
264  tempshift.SetShift(systConfig.first, (double) rand.Gaus());
265  }
266  else if (systConfig.second == kSystInteger){
267  tempshift.SetShift(systConfig.first, (double) rand.Integer(2));
268  }
269  } //end systconfigs
270  shifts.push_back(tempshift);
271  }//end universes
272  return shifts;
273  }
274 
275  //----------------------------------------------------------------------
276  std::unique_ptr<SystShifts> GaussianPriorSystShifts::Copy() const
277  {
278  return std::make_unique<GaussianPriorSystShifts>(*this);
279  }
280 
281  //----------------------------------------------------------------------
283  {
284  stan::math::var ret = 0;
285 
286  // Systematics are all expressed in terms of sigmas;
287  // the Gaussian prior means its log is quadratic...
288  for (const auto & systPair: fSystsStan)
289  ret -= 0.5 * util::sqr(systPair.second);
290 
291  return ret;
292  }
293 
294  //----------------------------------------------------------------------
296  {
297  stan::math::var ret = 0;
298 
299  // Systematics are all expressed in terms of sigmas
300  for (const auto & systPair: fSystsStan)
301  ret += exp(-0.5 * util::sqr(systPair.second));
302 
303  return ret;
304 
305  }
306 }
static int fgNextID
The next unused ID.
Definition: SystShifts.h:91
std::string ShortName() const
Brief description of component shifts, for printing to screen.
Definition: SystShifts.cxx:171
Cuts and Vars for the 2020 FD DiF Study.
Definition: CutFlow_header.h:5
virtual stan::math::var Prior() const
Definition: SystShifts.h:65
virtual std::unique_ptr< SystShifts > Copy() const
Definition: SystShifts.cxx:64
set< int >::iterator it
bool IsNominal() const
Definition: SystShifts.h:43
virtual const std::string & ShortName() const final
The name printed out to the screen.
Definition: ISyst.h:27
const Var weight
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
Proxy for caf::StandardRecord.
Definition: SRProxy.h:1969
static std::unique_ptr< SystShifts > LoadFrom(TDirectory *dir)
Definition: SystShifts.cxx:230
std::unordered_map< const ISyst *, stan::math::var > fSystsStan
Definition: SystShifts.h:86
OStream cerr
Definition: OStream.cxx:7
Float_t tmp
Definition: plot.C:36
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
std::vector< SystShifts > GetSystShiftsMultiverse(const std::vector< std::pair< const ISyst *, SystMode > > systConfigs, const int nUniverses, const int seed)
Definition: SystShifts.cxx:250
unsigned int seed
Definition: runWimpSim.h:102
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
std::string LatexName() const
Long description of component shifts, for plot labels.
Definition: SystShifts.cxx:185
T GetShift(const ISyst *syst) const
stan::math::var Prior() const override
Definition: SystShifts.cxx:295
std::unordered_map< const ISyst *, double > fSystsDbl
Definition: SystShifts.h:85
OStream cout
Definition: OStream.cxx:6
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
TDirectory * dir
Definition: macro.C:5
void ResetToNominal()
Definition: SystShifts.cxx:133
void SaveTo(TDirectory *dir) const
Definition: SystShifts.cxx:207
assert(nhit_max >=nhit_nbins)
virtual stan::math::var LogPrior() const
Definition: SystShifts.cxx:151
std::vector< const ISyst * > ActiveSysts() const
Definition: SystShifts.cxx:199
double Penalty() const
Penalty term for (frequentist) chi-squared fits.
Definition: SystShifts.cxx:142
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
stan::math::var LogPrior() const override
Definition: SystShifts.cxx:282
std::unique_ptr< SystShifts > Copy() const override
Definition: SystShifts.cxx:276
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:70
static constexpr Double_t sr
Definition: Munits.h:164
void Shift(caf::SRProxy *sr, double &weight) const
Definition: SystShifts.cxx:154