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 #include <math.h>
12 
13 #include "TDirectory.h"
14 #include "TH1.h"
15 #include "TMath.h"
16 #include "TObjString.h"
17 #include "TString.h"
18 #include "TRandom3.h"
19 
20 namespace ana
21 {
23 
24  // Reserve 0 for unshifted
25  int SystShifts::fgNextID = 1;
26 
27  //----------------------------------------------------------------------
29  : fID(fgNextID++)
30  {}
31 
32  //----------------------------------------------------------------------
33  SystShifts::SystShifts(const ISyst* syst, double shift)
34  : fID(fgNextID++)
35  {
36  fSystsDbl.emplace(syst, shift);
37  }
38 
39  //----------------------------------------------------------------------
41  : fID(fgNextID++)
42  {
43  fSystsStan.emplace(syst, shift);
44  // we're always going to maintain a "double" copy
45  // so that when the Stan cache gets invalidated we can still return something usable.
46  fSystsDbl.emplace(syst, util::GetValAs<double>(shift));
47  }
48 
49  //----------------------------------------------------------------------
50  SystShifts::SystShifts(const std::map<const ISyst*, double>& shifts)
51  : fID(fgNextID++)
52  {
53  for(auto it: shifts) fSystsDbl.emplace(it.first, it.second);
54  }
55 
56  //----------------------------------------------------------------------
57  SystShifts::SystShifts(const std::map<const ISyst*, stan::math::var>& shifts)
58  : fID(fgNextID++)
59  {
60  for(auto it: shifts)
61  {
62  fSystsStan.emplace(it.first, it.second);
63  fSystsDbl.emplace(it.first, util::GetValAs<double>(it.second));
64  }
65  }
66 
67  //----------------------------------------------------------------------
68  std::unique_ptr<SystShifts> SystShifts::Copy() const
69  {
70  return std::make_unique<SystShifts>(*this);
71  }
72 
73  //----------------------------------------------------------------------
74  void SystShifts::RemoveShift(const ISyst *syst)
75  {
76  fSystsStan.erase(syst);
77  fSystsDbl.erase(syst);
78  }
79 
80  //----------------------------------------------------------------------
81  void SystShifts::SetShift(const ISyst* syst, double shift, bool force)
82  {
83  fID = fgNextID++;
84 
85  // if this slot already exists in the Stan systs, and we're not setting the same value,
86  // some shenanigans are going on that we need to figure out. abort.
87  auto itStan = fSystsStan.find(syst);
88  if(itStan != fSystsStan.end() && itStan->second != shift)
89  {
90  std::cerr << "Error Syst '" << syst->ShortName() << " already has a Stan pull set (" << itStan->second << ") "
91  << "and you're trying to set a different double one (" << shift << ")." << std::endl;
92  std::cerr << "You almost certainly didn't mean to do that." << std::endl;
93  std::cerr << "Abort." << std::endl;
94  abort();
95  }
96 
97  fSystsDbl.erase(syst);
98  if(force || shift != 0.) fSystsDbl.emplace(syst, shift);
99  }
100 
101  //----------------------------------------------------------------------
102  void SystShifts::SetShift(const ISyst* syst, stan::math::var shift)
103  {
104  fSystsStan.erase(syst);
105  // note: _always_ put the syst in, even if the value is 0.
106  // autodiff relies on the calculation happening so as to get the gradient
107  fSystsStan.emplace(syst, shift);
108  SetShift(syst, util::GetValAs<double>(shift), true);
109  }
110 
111  //----------------------------------------------------------------------
112  template <>
113  double SystShifts::GetShift(const ISyst* syst) const
114  {
115  assert(syst);
116 
117  auto it = fSystsDbl.find(syst);
118  return (it == fSystsDbl.end()) ? 0 : it->second;
119  }
120 
121  //----------------------------------------------------------------------
122  template <>
124  {
125  assert(syst);
126 
127  auto it = fSystsStan.find(syst);
128  //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
129  if (it == fSystsStan.end())
130  {
131  auto itDbl = fSystsDbl.find(syst);
132  if (itDbl != fSystsDbl.end())
133  {
134  std::cout << "Warning: creating stan::math::var out of double value for syst '" << syst->ShortName() << "'." << std::endl;
135  std::cout << " If you see this repeatedly, it's likely a problem. (A few times during startup is probably harmless.)" << std::endl;
136  fSystsStan[syst] = stan::math::var(itDbl->second);
137  it = fSystsStan.find(syst);
138  }
139  }
140  return (it == fSystsStan.end()) ? 0 : it->second;
141  }
142 
143  //----------------------------------------------------------------------
145  {
146  fID = 0;
147 
148  fSystsDbl.clear();
149  fSystsStan.clear();
150  }
151 
152  //----------------------------------------------------------------------
153  double SystShifts::Penalty() const
154  {
155  double ret = 0;
156  // Systematics are all expressed in terms of sigmas
157  for(auto it: fSystsDbl) ret += it.second * it.second;
158  return ret;
159  }
160 
161  //----------------------------------------------------------------------
163 
164  //----------------------------------------------------------------------
166  double& weight) const
167  {
168  // always fSystsDbl here because this is only used in the event loop, not in fitting
169  // (so the autodiff'd version is not needed)
170  for(auto it: fSystsDbl){
171  double w = 1;
172  it.first->Shift(it.second, sr, w);
173  if(w < 0){
174  std::cout << "SystShifts: Got negative weight " << w
175  << " while shifting "
176  << it.first->ShortName()
177  << " to " << it.second << " sigma " << std::endl;
178  abort();
179  }
180  weight *= w;
181  }
182  }
183 
184  //----------------------------------------------------------------------
186  {
187  // always fSystsDbl here because this is only used in the event loop, not in fitting
188  // (so the autodiff'd version is not needed)
189  for(auto it: fSystsDbl) it.first->TruthShift(it.second, nu, weight);
190  }
191 
192  //----------------------------------------------------------------------
194  {
195  if(IsNominal()) return "nominal";
196 
198  for(auto it: fSystsDbl){
199  if(!ret.empty()) ret += ",";
200  ret += it.first->ShortName() + TString::Format("=%+g", it.second).Data();
201  }
202 
203  return ret;
204  }
205 
206  //----------------------------------------------------------------------
208  {
209  if(IsNominal()) return "Nominal";
210 
212  for(auto it: fSystsDbl){
213  if(!ret.empty()) ret += ", ";
214  ret += it.first->LatexName() + TString::Format(" = %+g", it.second).Data();
215  }
216 
217  return ret;
218  }
219 
220  //----------------------------------------------------------------------
221  std::vector<const ISyst*> SystShifts::ActiveSysts() const
222  {
223  std::vector<const ISyst*> ret;
224  for(auto it: fSystsDbl) ret.push_back(it.first);
225  return ret;
226  }
227 
228  //----------------------------------------------------------------------
229  void SystShifts::SaveTo(TDirectory* dir, const std::string& name) const
230  {
231  TDirectory* tmp = gDirectory;
232 
233  dir = dir->mkdir(name.c_str()); // switch to subdir
234  dir->cd();
235 
236  TObjString("SystShifts").Write("type");
237 
238  // Don't write any histogram for the nominal case
239  if(!fSystsDbl.empty()){
240  TH1D h("", "", fSystsDbl.size(), 0, fSystsDbl.size());
241  int ibin = 0;
242  for(auto it: fSystsDbl){
243  ++ibin;
244  h.GetXaxis()->SetBinLabel(ibin, it.first->ShortName().c_str());
245  h.SetBinContent(ibin, it.second);
246  }
247  h.Write("vals");
248  }
249 
250  dir->Write();
251  delete dir;
252 
253  tmp->cd();
254  }
255 
256  //----------------------------------------------------------------------
257  std::unique_ptr<SystShifts> SystShifts::LoadFrom(TDirectory* dir, const std::string& name)
258  {
259  dir = dir->GetDirectory(name.c_str()); // switch to subdir
260  assert(dir);
261 
262  TObjString* tag = (TObjString*)dir->Get("type");
263  assert(tag);
264  assert(tag->GetString() == "SystShifts");
265  delete tag;
266 
267  auto ret = std::make_unique<SystShifts>();
268 
269  TH1* h = (TH1*)dir->Get("vals");
270  if(h){ // no histogram means nominal
271  for(int i = 1; i <= h->GetNbinsX(); ++i){
272  ret->SetShift(Registry<ISyst>::ShortNameToPtr(h->GetXaxis()->GetBinLabel(i)),
273  h->GetBinContent(i));
274  }
275  }
276 
277  delete dir;
278 
279  return ret;
280  }
281 
282  //----------------------------------------------------------------------
283  std::vector <SystShifts> GetSystShiftsMultiverse(
284  const std::vector <std::pair <const ISyst*, SystMode> >& systConfigs,
285  const int nUniverses, const int seed)
286  {
287  std::vector <SystShifts> shifts = {};
288  int tempseed = seed;
289 
290  for (int univ=0; univ<nUniverses; ++univ){
291  SystShifts tempshift;
292 
293  for (auto systConfig:systConfigs) {
294  TRandom3 rand (tempseed++);
295  if(systConfig.second == kSystExclude) continue;
296  else if (systConfig.second == kSystGaussian) {
297  tempshift.SetShift(systConfig.first, (double) rand.Gaus());
298  }
299  else if (systConfig.second == kSystInteger){
300  tempshift.SetShift(systConfig.first, (double) rand.Integer(2));
301  }
302  } //end systconfigs
303  shifts.push_back(tempshift);
304  }//end universes
305  return shifts;
306  }
307 
308  //----------------------------------------------------------------------
309  std::unique_ptr<SystShifts> GaussianPriorSystShifts::Copy() const
310  {
311  return std::make_unique<GaussianPriorSystShifts>(*this);
312  }
313 
314  //----------------------------------------------------------------------
316  {
317  stan::math::var ret = 0;
318  // Systematics are all expressed in terms of sigmas
319  const double norm = log(2 * TMath::Pi()); // assumes sigma = 1, which is ok b/c we're working in units of sigma
320  for(auto it: fSystsStan) ret -= 0.5 * (norm + it.second * it.second);
321  return ret;
322  }
323 
324  //----------------------------------------------------------------------
326  {
327  return exp(LogPrior());
328  }
329 
330  //----------------------------------------------------------------------
331  std::unique_ptr<SystShifts> BoundedGaussianPriorSystShifts::Copy() const
332  {
333  return std::make_unique<BoundedGaussianPriorSystShifts>(*this);
334  }
335 
336  //----------------------------------------------------------------------
338  {
339  //Parabolic Prior shape, but with strict boundary at +/-3 sigma.
340  stan::math::var ret = 0;
341  const double norm = log(2 * TMath::Pi()); // assumes sigma = 1, which is ok b/c we're working in units of sigma
342  for(auto it: fSystsStan)
343  {
344  if (it.second < -3.0 || it.second > 3.0)
345  {
346  const double adjustment = 8094.083927575384; // difference to match the fns at 3sigma.
347  ret -= exp( (it.second * it.second) ) - adjustment;
348  }
349  else
350  {
351  ret -= 0.5 * (norm + it.second * it.second);
352  }
353  }
354  return ret;
355  }
356 
358  {
359  return exp(LogPrior());
360  }
361 } // namespace ana
const XML_Char * name
Definition: expat.h:151
static int fgNextID
The next unused ID.
Definition: SystShifts.h:96
std::string ShortName() const
Brief description of component shifts, for printing to screen.
Definition: SystShifts.cxx:193
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
virtual stan::math::var Prior() const
Definition: SystShifts.h:72
virtual std::unique_ptr< SystShifts > Copy() const
Definition: SystShifts.cxx:68
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
Proxy for caf::SRNeutrino.
Definition: SRProxy.h:510
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
std::unique_ptr< SystShifts > Copy() const override
Definition: SystShifts.cxx:331
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
std::unordered_map< const ISyst *, stan::math::var > fSystsStan
Definition: SystShifts.h:91
OStream cerr
Definition: OStream.cxx:7
Float_t tmp
Definition: plot.C:36
static SystShifts Nominal()
Definition: SystShifts.h:34
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
static std::unique_ptr< SystShifts > LoadFrom(TDirectory *dir, const std::string &name)
Definition: SystShifts.cxx:257
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:207
T GetShift(const ISyst *syst) const
caf::StandardRecord * sr
stan::math::var Prior() const override
Definition: SystShifts.cxx:325
std::vector< SystShifts > GetSystShiftsMultiverse(const std::vector< std::pair< const ISyst *, SystMode > > &systConfigs, const int nUniverses, const int seed)
Definition: SystShifts.cxx:283
stan::math::var Prior() const override
Definition: SystShifts.cxx:357
std::unordered_map< const ISyst *, double > fSystsDbl
Definition: SystShifts.h:90
const SystShifts kNoShift
Definition: SystShifts.cxx:22
OStream cout
Definition: OStream.cxx:6
Float_t norm
TDirectory * dir
Definition: macro.C:5
void ResetToNominal()
Definition: SystShifts.cxx:144
assert(nhit_max >=nhit_nbins)
virtual stan::math::var LogPrior() const
Definition: SystShifts.cxx:162
std::vector< const ISyst * > ActiveSysts() const
Definition: SystShifts.cxx:221
double Penalty() const
Penalty term for (frequentist) chi-squared fits.
Definition: SystShifts.cxx:153
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: SystShifts.cxx:229
stan::math::var LogPrior() const override
Definition: SystShifts.cxx:315
void RemoveShift(const ISyst *syst)
Definition: SystShifts.cxx:74
Float_t w
Definition: plot.C:20
std::unique_ptr< SystShifts > Copy() const override
Definition: SystShifts.cxx:309
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:81
void Shift(caf::SRProxy *sr, double &weight) const
Definition: SystShifts.cxx:165
stan::math::var LogPrior() const override
Definition: SystShifts.cxx:337
enum BeamMode string