MultiExperiment.cxx
Go to the documentation of this file.
2 
3 #include "CAFAna/Core/ISyst.h"
6 
7 #include "TDirectory.h"
8 #include "TObjString.h"
9 
10 #include <cassert>
11 #include <iostream>
12 
13 namespace ana
14 {
15  REGISTER_LOADFROM("MultiExperiment", IExperiment, MultiExperiment);
16 
17  //----------------------------------------------------------------------
19  const std::vector<std::pair<const ISyst *, const ISyst *>> &corrs) const
20  {
21  // Make a local copy we're going to rewrite into the terms this
22  // sub-experiment will accept.
23  SystShifts localShifts = shiftsIn;
24  for(auto it: corrs)
25  {
26  // We're mapping prim -> sec
27  const ISyst* prim = it.first;
28  const ISyst* sec = it.second;
29  if(shiftsIn.HasStan(prim) || shiftsIn.GetShift<double>(prim) != 0)
30  {
31  // sec can be unset, which means there's no representation needed
32  // of prim in the sub-experiment.
33  if(sec)
34  {
35  if (shiftsIn.HasStan(prim))
36  localShifts.SetShift(sec, shiftsIn.GetShift<stan::math::var>(prim));
37  else
38  localShifts.SetShift(sec, shiftsIn.GetShift<double>(prim));
39  }
40 
41  // We've either translated or discarded prim, so drop it here.
42  localShifts.RemoveShift(prim);
43  }
44  }
45 
46  return localShifts;
47  }
48 
49  //----------------------------------------------------------------------
51  const SystShifts& syst) const
52  {
53  double ret = 0;
54  for(unsigned int n = 0; n < fExpts.size(); ++n){
55  // Don't waste time fiddling with the systematics if for sure there's
56  // nothing to do.
57  if(fSystCorrelations[n].empty())
58  ret += fExpts[n]->ChiSq(osc, syst);
59  else
60  ret += fExpts[n]->ChiSq(osc, ApplySystCorrs(syst, fSystCorrelations[n]));
61  }
62  return ret;
63  }
64 
65  //----------------------------------------------------------------------
68  const std::vector<std::pair<const ISyst*, const ISyst*>>& corrs)
69  {
70  // Sanity-check the mapping
71  std::map<const ISyst*, const ISyst*> already;
72  for(auto it: corrs){
73  if(it.first == it.second)
74  {
75  std::cout << "MultiExperiment::SetSystCorrelations(): You're trying to map a syst to itself: " << std::endl;
76  auto name = it.first ? it.first->ShortName() : "nullptr";
77  std::cout << name << " --> " << name << std::endl;
78  abort();
79  }
80 
81  // Don't worry if second element is null pointer
82  if (!it.second) continue;
83  if(already.find(it.second) == already.end()){
84  already[it.second] = it.first;
85  }
86  else{
87  std::cout << "MultiExperiment::SetSystCorrelations(): Warning!\n"
88  << "In experiment " << idx << " both "
89  << already[it.second]->ShortName() << " and "
90  << it.first->ShortName()
91  << " are configured to map to " << it.second->ShortName()
92  << ". That's probably not what you want." << std::endl;
93  }
94  }
95 
96  // Apply it
97  fSystCorrelations[idx] = corrs;
98  }
99 
100  //----------------------------------------------------------------------
103  {
104  stan::math::var ret = 0.;
105  for (unsigned int idx = 0; idx < fExpts.size(); ++idx)
106  {
107  if(fSystCorrelations[idx].empty())
108  ret += fExpts[idx]->LogLikelihood(osc, syst);
109  else
110  ret += fExpts[idx]->LogLikelihood(osc, ApplySystCorrs(syst, fSystCorrelations[idx]));
111  }
112  return ret;
113  }
114 
115  //----------------------------------------------------------------------
116  void MultiExperiment::SaveTo(TDirectory* dir, const std::string& name) const
117  {
118  bool hasCorr = false;
119  for(auto it: fSystCorrelations) if(!it.empty()) hasCorr = true;
120 
121  if(hasCorr){
122  std::cerr << "Warning in MultiExperiment: systematic correlations are set and will not be serialized by this call to SaveTo(). You will have to re-set them once you load the experiment back in." << std::endl;
123  }
124 
125  TDirectory* tmp = dir;
126 
127  dir = dir->mkdir(name.c_str()); // switch to subdir
128  dir->cd();
129 
130  TObjString("MultiExperiment").Write("type");
131 
132  for(unsigned int i = 0; i < fExpts.size(); ++i){
133  fExpts[i]->SaveTo(dir, TString::Format("expt%d", i).Data());
134  }
135 
136  dir->Write();
137  delete dir;
138 
139  tmp->cd();
140  }
141 
142  //----------------------------------------------------------------------
143  std::unique_ptr<MultiExperiment> MultiExperiment::LoadFrom(TDirectory* dir, const std::string& name)
144  {
145  dir = dir->GetDirectory(name.c_str()); // switch to subdir
146  assert(dir);
147 
148  TObjString* ptag = (TObjString*)dir->Get("type");
149  assert(ptag);
150  assert(ptag->GetString() == "MultiExperiment");
151  delete ptag;
152 
153  std::vector<const IExperiment*> expts;
154 
155  for(int i = 0; ; ++i){
156  const std::string subname = TString::Format("expt%d", i).Data();
157  TDirectory* subdir = dir->GetDirectory(subname.c_str());
158  if(!subdir) break;
159  delete subdir;
160 
161  expts.push_back(ana::LoadFrom<IExperiment>(dir, name).release());
162  }
163 
164  assert(!expts.empty());
165 
166  delete dir;
167 
168  return std::unique_ptr<MultiExperiment>(new MultiExperiment(expts));
169  }
170 
171 }
const XML_Char * name
Definition: expat.h:151
std::vector< const IExperiment * > fExpts
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
set< int >::iterator it
subdir
Definition: cvnie.py:7
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
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
std::vector< std::vector< std::pair< const ISyst *, const ISyst * > > > fSystCorrelations
stan::math::var LogLikelihood(osc::IOscCalcAdjustableStan *osc, const SystShifts &syst) const override
Sum up log-likelihoods of sub-expts.
T GetShift(const ISyst *syst) const
virtual double ChiSq(osc::IOscCalcAdjustable *osc, const SystShifts &syst=SystShifts::Nominal()) const override
bool HasStan(const ISyst *s) const
Definition: SystShifts.h:63
Oscillation probability calculators.
Definition: Calcs.h:5
OStream cout
Definition: OStream.cxx:6
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
SystShifts ApplySystCorrs(const SystShifts &shiftsIn, const std::vector< std::pair< const ISyst *, const ISyst * >> &corrs) const
Get a new SystShifts, applying the specified correlations to it.
TDirectory * dir
Definition: macro.C:5
REGISTER_LOADFROM("BENDecomp", IDecomp, BENDecomp)
assert(nhit_max >=nhit_nbins)
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
virtual void SaveTo(TDirectory *dir, const std::string &name) const override
std::unique_ptr< IExperiment > LoadFrom< IExperiment >(TDirectory *dir, const std::string &label)
Definition: IExperiment.cxx:16
static std::unique_ptr< MultiExperiment > LoadFrom(TDirectory *dir, const std::string &name)
void RemoveShift(const ISyst *syst)
Definition: SystShifts.cxx:73
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:80
MultiExperiment(std::vector< const IExperiment * > expts={})
void SetSystCorrelations(int idx, const std::vector< std::pair< const ISyst *, const ISyst * >> &corrs)