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"
15 
17 
18 using namespace std;
20 
21 #include "OscLib/IOscCalc.h"
22 
23 namespace mcmc
24 {
25  // ------------------------------------------------------------------------------
26 
27  std::vector<std::unique_ptr<ana::IExperiment>>
28  BuildExperiments(const std::map<std::string, ana::Spectrum>& fakeData,
29  const std::vector<ana::predictions>& preds)
30  {
31  std::vector<std::unique_ptr<ana::IExperiment>> ret;
32  for (const auto & predBundle : preds)
33  {
34  if ( fakeData.find(predBundle.name) == fakeData.end() )
35  {
36  std::cerr << "Couldn't find fake data for prediction: " << predBundle.name << std::endl;
37  abort();
38  }
39 
40  if (predBundle.cos.first)
41  ret.emplace_back(std::make_unique<ana::SingleSampleExperiment>(predBundle.pred,
42  fakeData.at(predBundle.name),
43  *predBundle.cos.first,
44  predBundle.cos.second));
45  else
46  ret.emplace_back(std::make_unique<ana::SingleSampleExperiment>(predBundle.pred,
47  fakeData.at(predBundle.name)));
48  }
49 
50  return ret;
51  }
52 
53  // ------------------------------------------------------------------------------
54  std::unique_ptr<ana::IExperiment> BuildMultiExperiment(const std::vector<const ana::IExperiment*>& expts)
55  {
56  auto multiexpt = std::make_unique<ana::MultiExperiment>();
57  std::set<const ana::ISyst*> allSysts;
58  for (const auto & expt : expts)
59  {
60  auto thisExptSysts = mcmc::SupportedSysts(expt);
61  std::copy(thisExptSysts.begin(), thisExptSysts.end(), std::inserter(allSysts, allSysts.end()));
62  }
63 
64  for (std::size_t exptIdx = 0; exptIdx < expts.size(); exptIdx++)
65  {
66  auto expt = expts[exptIdx];
67  multiexpt->Add(expt);
68 
69  auto thisExptSysts = mcmc::SupportedSysts(expt);
70  std::vector<std::pair<const ana::ISyst*, const ana::ISyst*>> correlations;
71  for (const auto & syst : allSysts)
72  {
73  if (std::find(thisExptSysts.begin(), thisExptSysts.end(), syst) == thisExptSysts.end())
74  correlations.emplace_back(std::make_pair(syst, nullptr));
75  }
76  multiexpt->SetSystCorrelations(exptIdx, correlations);
77  }
78 
79  return multiexpt;
80  }
81 
82  // ------------------------------------------------------------------------------
83 
84  std::map<std::string, ana::Spectrum>
85  LoadFakeData(const std::string & fakeDataFilename, const std::string & path)
86  {
87  std::map<std::string, ana::Spectrum> ret;
88 
89  std::cout << "loading from fake data file: " << fakeDataFilename << std::endl;
90  TFile inF(fakeDataFilename.c_str());
91  if (inF.IsZombie())
92  abort();
93 
94  TDirectory * workDir = &inF;
95  if (!path.empty())
96  workDir = inF.GetDirectory(path.c_str());
97  if (!workDir)
98  {
99  std::cerr << "Couldn't open path '" << path << "' in file: '" << inF.GetName() << std::endl;
100  abort();
101  }
102 
103  for (const auto & k : *workDir->GetListOfKeys())
104  {
105  auto key = dynamic_cast<TKey*>(k);
106  auto subdir = dynamic_cast<TDirectory*>(key->ReadObj());
107  if (!subdir)
108  continue;
109 
110  if (!subdir->Get("type")
111  || dynamic_cast<TObjString*>(subdir->Get("type"))->GetString() != "Spectrum")
112  continue;
113 
114  std::cout << "Loading fake data spectrum: '" << key->GetName() << "'" << std::endl;
115  auto spec = ana::LoadFrom<ana::Spectrum>(workDir, key->GetName());
116  if (spec)
117  ret.emplace(std::piecewise_construct,
118  std::forward_as_tuple(key->GetName()),
119  std::forward_as_tuple(*spec.release()));
120  }
121 
122  return ret;
123  }
124 
125  // ------------------------------------------------------------------------------
126 
128  {
129  TFile inf(filename.c_str());
130  if (inf.IsZombie())
131  return {};
132 
133  std::unique_ptr<ana::Spectrum> fakeData;
134  if (inf.Get(config::FAKE_DATA_LABEL.c_str()))
135  fakeData = ana::LoadFrom<ana::Spectrum>(&inf, config::FAKE_DATA_LABEL);
136  std::unique_ptr<ana::MCMCSamples> warmup, samples;
137  if (inf.Get(config::WARMUP_LABEL.c_str()))
138  warmup = ana::LoadFrom<ana::MCMCSamples>(&inf, config::WARMUP_LABEL);
139  if (inf.Get(config::SAMPLES_LABEL.c_str()))
140  samples = ana::LoadFrom<ana::MCMCSamples>(&inf, config::SAMPLES_LABEL);
141  std::unique_ptr<ana::SystShifts> systTruePulls;
142  if (inf.Get(config::SYST_PULLS_LABEL.c_str()))
143  systTruePulls = ana::LoadFrom<ana::SystShifts>(&inf, config::SYST_PULLS_LABEL);
144  std::unique_ptr<osc::IOscCalcAdjustable> trueOscCalc;
145  if (inf.Get(config::TRUE_OSC_LABEL.c_str()))
146  trueOscCalc = ana::LoadFrom<osc::IOscCalcAdjustable>(&inf, config::TRUE_OSC_LABEL);
147 
148  return FileContents
149  {
150  std::move(fakeData),
151  std::move(warmup),
152  std::move(samples),
153  std::move(systTruePulls),
154  std::move(trueOscCalc)
155  };
156  }
157 
158  // ------------------------------------------------------------------------------
159  //previously LoadPreds()
160  std::vector<ana::predictions> LoadStock3FFDPreds(bool extrap)
161  {
162  if (extrap)
163  std::cerr << "MCMC3FShared.cxx::LoadStock3FFDPreds(): Warning: extrapolated predictions not implemented yet" << std::endl;
164  return ana::LoadPredictions(true /* load systs? */,
165  false, /* on grid? --> not actually used currently */
166  extrap ? "combo" : "noextrap", /* decomp desired. empty string means default: full |pt| extrap.
167  can indicate no-|pt| using "noPt".
168  combine this with "noextrap", "combo", or "prop"
169  to specify nue decomp if desired */
170  true, /* load nue predictions? */
171  true, /* load numu predictions? */
172  true, /* load FHC predictions? */
173  true, /* load RHC predictions? */
174  false //, /* running on NERSC? */
175 // extrap
176  );
177 
178  }
179  // ------------------------------------------------------------------------------
180  std::vector<ana::predictions> LoadNDPreds(const std::vector<std::pair<std::string, std::string>> & filenames,
181  const std::vector<std::string> & predObjNames)
182  {
183 // // for now, bc we are testing the AllNumu of the ND predictions file:
184 // std::vector<std::string> predObjNames{
185 // "pred_interp_AllNumu"
186 // };
187 
188  std::vector<ana::predictions> preds;
189 
190 
191  for (const auto & filenamePair : filenames)
192  {
193  const std::string sampleName = filenamePair.first;
194  // open fhc file
195  TFile f(filenamePair.second.c_str());
196 
197  // check file is valid
198  if (f.IsZombie())
199  {
200  std::cerr << "Prediction not found for: " << filenamePair.first << std::endl;
201  abort();
202  }
203 
204 
205  for (const auto & predName : predObjNames)
206  {
207  auto predInterp = ana::LoadFrom<ana::PredictionInterp>(&f, predName);
208  if (!predInterp)
209  {
210  // no prediction loaded from file
211  std::cerr << "File not not found. LoadFrom() failed to find file: " << &f << std::endl;
212  exit(1);
213  }
214  ana::predictions pred //temporary pred
215  {
216  sampleName + "_" + predName,
217  predInterp.release(),
218  {nullptr, 0},
219  ana::kAna2020FHCFullDetEquivPOT,// <-- todo: fix
220  ana::kAna2020FHCLivetime // <-- todo: fix this
221  };
222  preds.emplace_back(std::move(pred));
223  }
224  }
225  return preds;
226  }
227 
228 // ------------------------------------------------------------------------------
229  std::vector<std::unique_ptr<ana::ISyst>> LoadSystsCreateDummyMECSysts()
230  {
231  //load all xsec systs, but do not save output
233  //create dummy systs for now on the MEC knobs
234  std::vector<std::string> knobsToFake
235  {
236  "MECDoubleGaussEnhSystSigmaQ0_1",
237  "MECDoubleGaussEnhSystCorr_1",
238  "MECDoubleGaussEnhSystMeanQ0_1",
239  "MECDoubleGaussEnhSystMeanQ3_2",
240  "MECDoubleGaussEnhSystSigmaQ3_2",
241  "MECDoubleGaussEnhSystNorm_1",
242  "MECDoubleGaussEnhSystNorm_2",
243  "MECDoubleGaussEnhSystMeanQ3_1",
244  "MECDoubleGaussEnhSystCorr_2",
245  "MECDoubleGaussEnhSystMeanQ0_2",
246  "MECDoubleGaussEnhSystSigmaQ0_2",
247  "MECDoubleGaussEnhSystSigmaQ3_1",
248  "MECDoubleGaussEnhSystBaseline",
249  };
250  std::vector<std::unique_ptr<ana::ISyst>> ret;
251  for (std::size_t idx = 0; idx < knobsToFake.size(); idx++)
252  ret.emplace_back(std::make_unique<ana::DummyAnaSyst>(knobsToFake[idx], "Dummy MEC syst #" + std::to_string(idx+1)));
253  std::cout << "2020 Systs loaded." << std::endl;
254  std::cout << "MEC knobs are now DummyAnaSyst." << std::endl;
255  return ret;
256  }
257 
258  // ------------------------------------------------------------------------------
260  {
261  calc.SetL(810);
262  calc.SetRho(2.75);
263  calc.SetDmsq21(7.6e-5);
264  calc.SetDmsq32(2.5e-3);
265  calc.SetTh12(asin(sqrt(.87))/2);
266  calc.SetTh13(asin(sqrt(.10))/2);
267  calc.SetTh23(TMath::Pi() / 4);
268  calc.SetdCP(0);
269  }
270 
271  // ------------------------------------------------------------------------------
272 
274  const ana::Spectrum * fakeData,
275  const ana::MCMCSamples * warmup,
276  const ana::MCMCSamples * samples,
277  const ana::SystShifts * trueSystPulls,
278  const osc::IOscCalcAdjustable * trueOscParams)
279  {
280  TFile outF(fileName.c_str(), "recreate");
281  if (outF.IsZombie())
282  abort();
283 
284  if (fakeData)
285  fakeData->SaveTo(&outF, config::FAKE_DATA_LABEL);
286  if (warmup)
287  warmup->SaveTo(&outF, config::WARMUP_LABEL);
288  if (samples)
289  samples->SaveTo(&outF, config::SAMPLES_LABEL);
290  if (trueSystPulls)
291  trueSystPulls->SaveTo(&outF, config::SYST_PULLS_LABEL);
292  if (trueOscParams)
293  ana::SaveTo(*dynamic_cast<const osc::IOscCalc*>(trueOscParams), &outF, config::TRUE_OSC_LABEL);
294 
295  outF.Close();
296  }
297 
298  // ------------------------------------------------------------------------------
299  std::set<const ana::ISyst*> SupportedSysts(const ana::IExperiment* expt)
300  {
301  std::set<const ana::ISyst*> ret;
302  const ana::IPrediction * pred = nullptr;
303 
304  if (auto multiexpt = dynamic_cast<const ana::MultiExperiment*>(expt))
305  {
306  for (const auto &e : multiexpt->GetExperiments())
307  {
308  std::set<const ana::ISyst *> thisExptSysts = SupportedSysts(e);
309  std::copy(thisExptSysts.begin(), thisExptSysts.end(), std::inserter(ret, ret.end()));
310  }
311  }
312  if (auto ssexpt = dynamic_cast<const ana::SingleSampleExperiment*>(expt))
313  pred = ssexpt->GetPrediction();
314 
315  if (auto predAddRock = dynamic_cast<const ana::PredictionAddRock*>(pred))
316  {
317  ret.insert(&ana::kRockScaleSyst);
318  pred = predAddRock->GetFiducialPred();
319  }
320 
321  if (auto predInterp = dynamic_cast<const ana::PredictionInterp*>(pred))
322  {
323  auto interpSyst = predInterp->GetAllSysts();
324  std::copy(interpSyst.begin(), interpSyst.end(), std::inserter(ret, ret.end()));
325  }
326 
327  return ret;
328  }
329 
330 
331 }
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
std::vector< ana::predictions > LoadStock3FFDPreds(bool extrap)
Just wrap up all the args somewhere centralized for the moment.
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
std::vector< std::unique_ptr< ana::ISyst > > LoadSystsCreateDummyMECSysts()
osc::OscCalcDumb calc
virtual void SetDmsq32(const T &dmsq32)=0
const double kAna2020FHCFullDetEquivPOT
Definition: Exposures.h:234
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:39
std::vector< const ISyst * > getAllXsecSysts_2020()
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 double kAna2020FHCLivetime
Definition: Exposures.h:237
const std::string FAKE_DATA_LABEL
Definition: MCMC3FShared.h:45
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:534
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.
exit(0)
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:229
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 > LoadNDPreds(const std::vector< std::pair< std::string, std::string >> &filenames, const std::vector< std::string > &predObjNames)
std::string to_string(ModuleType const mt)
Definition: ModuleType.h:34
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