MCMC3FShared.cxx
Go to the documentation of this file.
2 
3 #include <iostream>
4 
5 #include "TKey.h"
6 #include "TMath.h"
7 
8 #include "CAFAna/Core/Spectrum.h"
13 #include "CAFAna/Fit/MCMCSamples.h"
14 
16 
17 using namespace std;
19 
20 #include "OscLib/IOscCalc.h"
21 
22 namespace mcmc
23 {
24  // ------------------------------------------------------------------------------
25 
26  std::vector<std::unique_ptr<ana::IExperiment>>
27  BuildExperiments(const std::map<std::string, ana::Spectrum>& fakeData,
28  const std::vector<ana::predictions>& preds)
29  {
30  std::vector<std::unique_ptr<ana::IExperiment>> ret;
31  for (const auto & predBundle : preds)
32  {
33  if ( fakeData.find(predBundle.name) == fakeData.end() )
34  {
35  std::cerr << "Couldn't find fake data for prediction: " << predBundle.name << std::endl;
36  abort();
37  }
38 
39  if (predBundle.cos.first)
40  ret.emplace_back(std::make_unique<ana::SingleSampleExperiment>(predBundle.pred,
41  fakeData.at(predBundle.name),
42  *predBundle.cos.first,
43  predBundle.cos.second));
44  else
45  ret.emplace_back(std::make_unique<ana::SingleSampleExperiment>(predBundle.pred,
46  fakeData.at(predBundle.name)));
47  }
48 
49  return ret;
50  }
51 
52  // ------------------------------------------------------------------------------
53  std::unique_ptr<ana::IExperiment> BuildMultiExperiment(const std::vector<const ana::IExperiment*>& expts)
54  {
55  auto multiexpt = std::make_unique<ana::MultiExperiment>();
56  std::set<const ana::ISyst*> allSysts;
57  for (const auto & expt : expts)
58  {
59  auto thisExptSysts = mcmc::SupportedSysts(expt);
60  std::copy(thisExptSysts.begin(), thisExptSysts.end(), std::inserter(allSysts, allSysts.end()));
61  }
62 
63  for (std::size_t exptIdx = 0; exptIdx < expts.size(); exptIdx++)
64  {
65  auto expt = expts[exptIdx];
66  multiexpt->Add(expt);
67 
68  auto thisExptSysts = mcmc::SupportedSysts(expt);
69  std::vector<std::pair<const ana::ISyst*, const ana::ISyst*>> correlations;
70  for (const auto & syst : allSysts)
71  {
72  if (std::find(thisExptSysts.begin(), thisExptSysts.end(), syst) == thisExptSysts.end())
73  correlations.emplace_back(std::make_pair(syst, nullptr));
74  }
75  multiexpt->SetSystCorrelations(exptIdx, correlations);
76  }
77 
78  return multiexpt;
79  }
80 
81  // ------------------------------------------------------------------------------
82 
83  std::map<std::string, ana::Spectrum>
84  LoadFakeData(const std::string & fakeDataFilename, const std::string & path)
85  {
86  std::map<std::string, ana::Spectrum> ret;
87 
88  std::cout << "loading from fake data file: " << fakeDataFilename << std::endl;
89  TFile inF(fakeDataFilename.c_str());
90  if (inF.IsZombie())
91  abort();
92 
93  TDirectory * workDir = &inF;
94  if (!path.empty())
95  workDir = inF.GetDirectory(path.c_str());
96  if (!workDir)
97  {
98  std::cerr << "Couldn't open path '" << path << "' in file: '" << inF.GetName() << std::endl;
99  abort();
100  }
101 
102  for (const auto & k : *workDir->GetListOfKeys())
103  {
104  auto key = dynamic_cast<TKey*>(k);
105  auto subdir = dynamic_cast<TDirectory*>(key->ReadObj());
106  if (!subdir)
107  continue;
108 
109  if (!subdir->Get("type")
110  || dynamic_cast<TObjString*>(subdir->Get("type"))->GetString() != "Spectrum")
111  continue;
112 
113  std::cout << "Loading fake data spectrum: '" << key->GetName() << "'" << std::endl;
114  auto spec = ana::LoadFrom<ana::Spectrum>(workDir, key->GetName());
115  if (spec)
116  ret.emplace(std::piecewise_construct,
117  std::forward_as_tuple(key->GetName()),
118  std::forward_as_tuple(*spec.release()));
119  }
120 
121  return ret;
122  }
123 
124  // ------------------------------------------------------------------------------
125 
127  {
128  TFile inf(filename.c_str());
129  if (inf.IsZombie())
130  return {};
131 
132  std::unique_ptr<ana::Spectrum> fakeData;
133  if (inf.Get(config::FAKE_DATA_LABEL.c_str()))
134  fakeData = ana::LoadFrom<ana::Spectrum>(&inf, config::FAKE_DATA_LABEL);
135  std::unique_ptr<ana::MCMCSamples> warmup, samples;
136  if (inf.Get(config::WARMUP_LABEL.c_str()))
137  warmup = ana::LoadFrom<ana::MCMCSamples>(&inf, config::WARMUP_LABEL);
138  if (inf.Get(config::SAMPLES_LABEL.c_str()))
139  samples = ana::LoadFrom<ana::MCMCSamples>(&inf, config::SAMPLES_LABEL);
140  std::unique_ptr<ana::SystShifts> systTruePulls;
141  if (inf.Get(config::SYST_PULLS_LABEL.c_str()))
142  systTruePulls = ana::LoadFrom<ana::SystShifts>(&inf, config::SYST_PULLS_LABEL);
143  std::unique_ptr<osc::IOscCalcAdjustable> trueOscCalc;
144  if (inf.Get(config::TRUE_OSC_LABEL.c_str()))
145  trueOscCalc = ana::LoadFrom<osc::IOscCalcAdjustable>(&inf, config::TRUE_OSC_LABEL);
146 
147  return FileContents
148  {
149  std::move(fakeData),
150  std::move(warmup),
151  std::move(samples),
152  std::move(systTruePulls),
153  std::move(trueOscCalc)
154  };
155  }
156 
157  // ------------------------------------------------------------------------------
158  std::vector<ana::predictions> LoadPreds(bool extrap)
159  {
160  if (extrap)
161  std::cerr << "MCMC3FShared.cxx::LoadPreds(): Warning: extrapolated predictions not implemented yet" << std::endl;
162  return ana::LoadPredictions(true /* load systs? */,
163  false, /* on grid? --> not actually used currently */
164  extrap ? "combo" : "noextrap", /* decomp desired. empty string means default: full |pt| extrap.
165  can indicate no-|pt| using "noPt".
166  combine this with "noextrap", "combo", or "prop"
167  to specify nue decomp if desired */
168  true, /* load nue predictions? */
169  true, /* load numu predictions? */
170  true, /* load FHC predictions? */
171  true, /* load RHC predictions? */
172  false //, /* running on NERSC? */
173 // extrap
174  );
175 
176  }
177 
178  // ------------------------------------------------------------------------------
179 
181  {
182  calc.SetL(810);
183  calc.SetRho(2.75);
184  calc.SetDmsq21(7.6e-5);
185  calc.SetDmsq32(2.5e-3);
186  calc.SetTh12(asin(sqrt(.87))/2);
187  calc.SetTh13(asin(sqrt(.10))/2);
188  calc.SetTh23(TMath::Pi() / 4);
189  calc.SetdCP(0);
190  }
191 
192  // ------------------------------------------------------------------------------
193 
195  const ana::Spectrum * fakeData,
196  const ana::MCMCSamples * warmup,
197  const ana::MCMCSamples * samples,
198  const ana::SystShifts * trueSystPulls,
199  const osc::IOscCalcAdjustable * trueOscParams)
200  {
201  TFile outF(fileName.c_str(), "recreate");
202  if (outF.IsZombie())
203  abort();
204 
205  if (fakeData)
206  fakeData->SaveTo(&outF, config::FAKE_DATA_LABEL);
207  if (warmup)
208  warmup->SaveTo(&outF, config::WARMUP_LABEL);
209  if (samples)
210  samples->SaveTo(&outF, config::SAMPLES_LABEL);
211  if (trueSystPulls)
212  trueSystPulls->SaveTo(&outF, config::SYST_PULLS_LABEL);
213  if (trueOscParams)
214  ana::SaveTo(*dynamic_cast<const osc::IOscCalc*>(trueOscParams), &outF, config::TRUE_OSC_LABEL);
215 
216  outF.Close();
217  }
218 
219  // ------------------------------------------------------------------------------
220  std::set<const ana::ISyst*> SupportedSysts(const ana::IExperiment* expt)
221  {
222  std::set<const ana::ISyst*> ret;
223  const ana::IPrediction * pred = nullptr;
224 
225  if (auto multiexpt = dynamic_cast<const ana::MultiExperiment*>(expt))
226  {
227  for (const auto &e : multiexpt->GetExperiments())
228  {
229  std::set<const ana::ISyst *> thisExptSysts = SupportedSysts(e);
230  std::copy(thisExptSysts.begin(), thisExptSysts.end(), std::inserter(ret, ret.end()));
231  }
232  }
233  if (auto ssexpt = dynamic_cast<const ana::SingleSampleExperiment*>(expt))
234  pred = ssexpt->GetPrediction();
235 
236  if (auto predAddRock = dynamic_cast<const ana::PredictionAddRock*>(pred))
237  {
238  ret.insert(&ana::kRockScaleSyst);
239  pred = predAddRock->GetFiducialPred();
240  }
241 
242  if (auto predInterp = dynamic_cast<const ana::PredictionInterp*>(pred))
243  {
244  auto interpSyst = predInterp->GetAllSysts();
245  std::copy(interpSyst.begin(), interpSyst.end(), std::inserter(ret, ret.end()));
246  }
247 
248  return ret;
249  }
250 
251 
252 }
virtual void SetL(double L)=0
fileName
Definition: plotROC.py:78
TFile * inf
Definition: drawXsec.C:1
void SaveTo(const osc::IOscCalc &x, TDirectory *dir, const std::string &name)
std::unique_ptr< ana::IExperiment > BuildMultiExperiment(const std::vector< const ana::IExperiment * > &expts)
virtual void SetDmsq21(const T &dmsq21)=0
subdir
Definition: cvnie.py:7
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
void SaveToFile(const std::string &fileName, const ana::Spectrum *fakeData, const ana::MCMCSamples *warmup, const ana::MCMCSamples *samples, const ana::SystShifts *trueSystPulls, const osc::IOscCalcAdjustable *trueOscParams)
T sqrt(T number)
Definition: d0nt_math.hpp:156
virtual void SetTh13(const T &th13)=0
OStream cerr
Definition: OStream.cxx:7
std::vector< predictions > LoadPredictions(bool corrSysts=false, bool runOnGrid=false, std::string decomp="", bool nue=true, bool numu=true, bool fhc=true, bool rhc=true, bool NERSC=false)
string filename
Definition: shutoffs.py:106
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
osc::OscCalcDumb calc
virtual void SetDmsq32(const T &dmsq32)=0
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
void SaveTo(TDirectory *dir, const std::string &name) const
Save this guy to a file so we don&#39;t have to rerun the MCMC.
expt
Definition: demo5.py:34
std::set< const ana::ISyst * > SupportedSysts(const ana::IExperiment *expt)
Which systs does this experiment support?
std::vector< std::unique_ptr< ana::IExperiment > > BuildExperiments(const std::map< std::string, ana::Spectrum > &fakeData, const std::vector< ana::predictions > &preds)
const std::string FAKE_DATA_LABEL
Definition: MCMC3FShared.h:45
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:506
OStream cout
Definition: OStream.cxx:6
const std::string path
Definition: plot_BEN.C:43
const std::string SYST_PULLS_LABEL
Definition: MCMC3FShared.h:47
virtual void SetRho(double rho)=0
string GetString(xmlDocPtr xml_doc, string node_path)
const std::string SAMPLES_LABEL
Definition: MCMC3FShared.h:46
virtual void SetTh23(const T &th23)=0
General interface to any calculator that lets you set the parameters.
Base class defining interface for experiments.
Definition: IExperiment.h:14
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
std::map< std::string, ana::Spectrum > LoadFakeData(const std::string &fakeDataFilename, const std::string &path)
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: SystShifts.cxx:228
const std::string TRUE_OSC_LABEL
Definition: MCMC3FShared.h:48
void ResetCalculator(osc::IOscCalcAdjustable &calc)
Float_t e
Definition: plot.C:35
std::vector< ana::predictions > LoadPreds(bool extrap)
Just wrap up all the args somewhere centralized for the moment.
FileContents LoadFromFile(const std::string &filename)
const std::string WARMUP_LABEL
Definition: MCMC3FShared.h:49
const DummyRockScaleSyst kRockScaleSyst
virtual void SetTh12(const T &th12)=0
virtual void SetdCP(const T &dCP)=0
T asin(T number)
Definition: d0nt_math.hpp:60
enum BeamMode string