ThrowFakeData.C
Go to the documentation of this file.
1 /*
2  * ThrowFakeData.C:
3  * Generate fake data for use in MCMC studies.
4  *
5  * Original author: J. Wolcott <jwolcott@fnal.gov>
6  * Date: September 2020
7  */
8 
9 #include "TRandom3.h"
10 
11 #include "OscLib/OscCalcPMNS.h"
12 #include "CAFAna/3flavor/Ana2020/joint_fit_2020_loader_tools.h"
13 
14 namespace mcmc_ana
15 {
16  // ---------------------------------------------
18  {
19  calc.SetL(810);
20  calc.SetRho(2.75);
21  calc.SetDmsq21(7.6e-5);
22  calc.SetDmsq32(2.5e-3);
23  calc.SetTh12(asin(sqrt(.87))/2);
24  calc.SetTh13(asin(sqrt(.10))/2);
25  calc.SetTh23(41 * TMath::DegToRad());
26  calc.SetdCP(5*TMath::Pi() / 4);
27  }
28 
29  // ---------------------------------------------
30  std::vector<const ana::ISyst*> GetSystList(const IPrediction* pred)
31  {
32  std::vector<const ana::ISyst*> systs;
33  if (auto predAddRock = dynamic_cast<const ana::PredictionAddRock*>(pred))
34  {
35  systs.emplace_back(&ana::kRockScaleSyst);
36  pred = predAddRock->GetFiducialPred();
37  }
38  if (auto predInterp = dynamic_cast<const ana::PredictionInterp*>(pred))
39  {
40  for (const auto & syst : predInterp->GetAllSysts())
41  systs.emplace_back(syst);
42  }
43  return systs;
44  }
45 
46 }
47 
48 void ThrowFakeData(const std::string & outFileName, bool throwSysts=true)
49 {
50  TFile outFile(outFileName.c_str(), "recreate");
51  if (outFile.IsZombie())
52  abort();
53 
56  ana::SaveTo(static_cast<osc::IOscCalc&>(calc), &outFile, "calcTruth"); // cast so the calculator-specific version gets called rather than the templated one
57 
58  // gimme 3-flavor predictions
59  std::vector<ana::predictions> preds = ana::LoadPredictions(true /* load systs? */,
60  false, /* on grid? --> not actually used currently */
61  "", /* decomp desired. empty string means default: full |pt| extrap.
62  can indicate no-|pt| using "noPt".
63  combine this with "noextrap", "combo", or "prop"
64  to specify nue decomp if desired */
65  true, /* load nue predictions? */
66  true, /* load numu predictions? */
67  true, /* load FHC predictions? */
68  true, /* load RHC predictions? */
69  false /* running on NERSC? */
70  );
71 
72  // construct a list of all of the possible systs with the pulls
73  std::unordered_set<const ana::ISyst*> allSysts;
74  for (const auto & predBundle : preds)
75  {
76  auto systs = mcmc_ana::GetSystList(predBundle.pred);
77  if (systs.empty())
78  std::cerr << "WARNING: prediction " << predBundle.name << " reports no systs..." << std::endl;
79  for (const auto & syst : systs)
80  allSysts.insert(syst);
81  }
82 
83  // choose pulls at random for all the possible systs.
84  TRandom3 rnd;
85  rnd.SetSeed(20200904);
86  ana::SystShifts shifts;
87  if (throwSysts)
88  {
89  for (const auto &syst : allSysts)
90  shifts.SetShift(syst, rnd.Gaus());
91  ana::SaveTo(shifts, &outFile, "systTruth");
92  }
93 
94  auto spectraDir = outFile.mkdir("spectra");
95  // now construct fake data spectra with the relevant systs
96  for (const auto predBundle : preds)
97  {
98  ana::SystShifts shiftsThisPred;
99  for (const auto & syst : mcmc_ana::GetSystList(predBundle.pred))
100  shiftsThisPred.SetShift(syst, shifts.GetShift(syst));
101  auto spec = predBundle.pred->PredictSyst(&calc, shiftsThisPred);
102  auto fakeData = spec.FakeData(predBundle.pot);
103  fakeData.SaveTo(spectraDir, predBundle.name.c_str());
104  }
105 
106  outFile.Close();
107 }
virtual void SetL(double L)=0
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
void SaveTo(const osc::IOscCalc &x, TDirectory *dir, const std::string &name)
virtual void SetDmsq21(const T &dmsq21)=0
Adapt the PMNS calculator to standard interface.
Definition: StanTypedefs.h:28
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
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)
osc::OscCalcDumb calc
virtual void SetDmsq32(const T &dmsq32)=0
void ResetCalculator(osc::IOscCalcAdjustable &calc)
Definition: ThrowFakeData.C:17
TFile * outFile
Definition: PlotXSec.C:135
T GetShift(const ISyst *syst) const
virtual void SetRho(double rho)=0
virtual void SetTh23(const T &th23)=0
std::vector< const ana::ISyst * > GetSystList(const IPrediction *pred)
Definition: ThrowFakeData.C:30
Float_t e
Definition: plot.C:35
void ThrowFakeData(const std::string &outFileName, bool throwSysts=true)
Definition: ThrowFakeData.C:48
#define for
Definition: msvc_pragmas.h:3
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:80
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