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