generate_fake_data.C
Go to the documentation of this file.
1 /// \file generate_fake_data.C
2 /// \brief Produce fake data for studying how different our 2019 results might look
3 /// given the 2018 data whose fit results we already know.
4 /// \author J. Wolcott <jwolcott@fnal.gov>
5 /// \date Feb. 7, 2019
6 
9 #include "CAFAna/Core/Spectrum.h"
10 #include "CAFAna/Core/Loaders.h"
12 
13 #include "fake_data.h"
14 
15 using namespace ana;
16 
17 namespace ana2019
18 {
19  /// Store the 2019 objects in tiered storage
20  template <typename StorageType>
22  {
23  public:
24  TieredAnaObjStorage(const std::function<const StorageType*(const std::string & beamStr)>& nueGetter,
25  const std::function<std::vector<const StorageType*>(const std::string & beamStr)>& numuGetter)
26  {
27  for (auto analysis : {kNueAna2018, kNumuAna2018})
28  {
29  for (auto beam : {kFHC, kRHC})
30  {
31  std::string beamStr = beam == kFHC ? "fhc" : "rhc";
32  if (analysis == kNueAna2018)
33  fStoredObjs[analysis][beam].emplace_back(std::unique_ptr<const StorageType>(nueGetter(beamStr)));
34  else if (analysis == kNumuAna2018)
35  {
36  auto objs = numuGetter(beamStr);
37  for (auto & obj : objs)
38  fStoredObjs[analysis][beam].emplace_back(std::unique_ptr<const StorageType>(obj));
39  }
40  }
41  }
42  }
43 
44  /// Note: quartileIfNumu counts from 1, not 0
45  const StorageType & Get(EAnaType2018 analysis, BeamType2018 beam, int quartileIfNumu=-1) const
46  {
47  assert (beam != kBoth);
48  assert (analysis != kJointAna2018);
49 
50  int idx = 0;
51  if (analysis == kNumuAna2018)
52  {
53  assert (quartileIfNumu >= 1 && quartileIfNumu <= 4);
54  idx = quartileIfNumu - 1;
55  }
56 
57  return *fStoredObjs.at(analysis).at(beam)[idx];
58  }
59 
60 
61  private:
62  std::map<EAnaType2018, std::map<BeamType2018, std::vector<std::unique_ptr<const StorageType>>>> fStoredObjs;
63  };
64 
66  [](const std::string & beamStr) { return GetNueData2018(beamStr); },
67  [](const std::string & beamStr)
68  {
69  auto vec = GetNumuData2018(4, beamStr);
70  std::vector<const Spectrum*> ret;
71  ret.insert(ret.end(), vec.begin(), vec.end());
72  return ret;
73  }
74  );
75 
77  // first bool is "use systs?"
78  [](const std::string & beamStr) { return GetNuePrediction2018(beamStr == "fhc" ? "combo" : "prop", DefaultOscCalc(), false, beamStr, false); },
79  [](const std::string & beamStr) { return GetNumuPredictions2018(4, false, beamStr); }
80  );
81 
82  // the GetNu*Cosmics() methods return pair<TH1, double> instead of Spectrum, so just do it manually
84  [](const std::string & beamStr)
85  {
87  std::string anaDir = "/nova/ana/nu_e_ana/Ana2018/";
88  if(beamStr=="rhc") name = "cosmic_spect_rhc";
89  else name="cosmic_spect_fhc";
90  auto spec = LoadFromFile<Spectrum>(anaDir+"/Predictions/cosmic/cosmic_prediction_real_data.root", name ).release();
91  return spec;
92  },
93 
94  [](const std::string & beamStr)
95  {
96  std::string dir = "/nova/ana/nu_mu_ana/Ana2018/Cosmics/";
97  std::string filename = (beamStr=="rhc") ? "cosmics_rhc__numu2018.root" : "cosmics_fhc__numu2018.root";
98  auto fcosm = TFile::Open((dir+filename).c_str());
99 
100  if(fcosm->IsZombie()) {std::cerr<< "bad cosmics\n"; exit(1);}
101  std::vector <const Spectrum*> numu_cosmics;
102 
103  for (int i = 1; i <=4; ++i)
104  numu_cosmics.emplace_back(new Spectrum(dynamic_cast<TH1*>(fcosm->Get(Form("cosmics_q%d", i))),
105  0, (beamStr == "rhc") ? kAna2018RHCLivetime : kAna2018FHCLivetime));
106 
107  return numu_cosmics;
108  }
109  );
110 
111  namespace fakedata
112  {
113 
114  /// Get a mock FD data spectrum based on the observed 2018 data spectra
115  /// plus predictions for the 2019 additional exposure based on the 2018 best fit
116  std::map<std::string, std::unique_ptr<Spectrum>>
118  EAnaType2018 analysis,
120  int quartileIfNumu=-1,
121  bool asimov=false)
122  {
123  std::map<std::string, std::unique_ptr<Spectrum>> ret;
124  ret.emplace("2018_data", std::make_unique<Spectrum>(kDataSpectra.Get(analysis, beam, quartileIfNumu)));
125 
126  // always start with the 2018 data
127  ret.emplace("2019_totalPred", std::make_unique<Spectrum>(*ret.at("2018_data")));
128 
129  // FHC will be unchanged for 2019
130  if (beam == kFHC)
131  return ret;
132 
133  std::cout << "analysis = " << analysis << ", beam = " << beam;
134  if (analysis == kNumuAna2018)
135  std::cout << ", quartile = " << quartileIfNumu;
136  std::cout << std::endl;
137 
138  // start with the data spectrum
139  auto & sumPred = ret.at("2019_totalPred");
140 
141  std::cout << " original data spectrum = " << ret.at("2018_data")->Integral(ret.at("2018_data")->POT(), nullptr, kPOT) << " events"
142  << " (" << kAna2018RHCPOT << " POT, " << kAna2018RHCLivetime << "s live)" << std::endl;
143 
144  // we are going to eventually have a spectrum at the new exposure.
145  // (anything we make below needs to have this same exposure, or it will wind up scaled to it...)
146  double scaleFactor = FULL_RHC_EXPOSURE / kAna2018RHCPOT;
147  sumPred->OverridePOT(FULL_RHC_EXPOSURE);
148 
149  // this cosmic spectrum was created from kAna2018RHCLivetime.
150  Spectrum cosmicsPred(kCosmicSpectra.Get(analysis, beam, quartileIfNumu));
151  cosmicsPred.OverridePOT(kAna2018RHCPOT * cosmicsPred.Livetime()/kAna2018RHCLivetime); // have to use POT for MockData()
152  std::cout << " (cosmics for 2018 RHC livetime: "
153  << cosmicsPred.Integral(kAna2018RHCLivetime, 0, kLivetime) << ")" << std::endl;
154 
155  double newPOT = (scaleFactor - 1) * kAna2018RHCPOT;
156  ret.emplace("2019_cosmicsPred", std::make_unique<Spectrum>(asimov ? cosmicsPred.FakeData(newPOT)
157  : cosmicsPred.MockData(newPOT)));
158  auto & cosmics = ret.at("2019_cosmicsPred");
159  // if its internal POT is different than that of 'ret'
160  // when we go to add it, this spectrum will be scaled so that
161  // its POT matches 'ret' before adding.
162  // since we've now produced the prediction at the correct exposure,
163  // we don't want it being rescaled, so just make its POT match
164  // and zero off its livetime so that its livetime isn't used at all.
165  cosmics->OverridePOT(FULL_RHC_EXPOSURE);
166  cosmics->OverrideLivetime(0);
167  *sumPred += *cosmics;
168 
169  std::cout << " cosmics = " << cosmics->Integral(cosmics->POT()) << " events"
170  << " (" << (scaleFactor - 1) * kAna2018RHCLivetime << "s live)" << std::endl;
171 
172  // finally, the beam prediction for the new exposure.
173  // this has to be overridden to the total exposure also
174  // so that it is added to 'ret' without rescaling
175  const IPrediction & pred = kPredictions.Get(analysis, beam, quartileIfNumu);
176  double newDataPOT = (scaleFactor-1)/scaleFactor * FULL_RHC_EXPOSURE;
177  if (asimov)
178  ret.emplace("2019_beamPred", std::make_unique<Spectrum>(pred.Predict(calc).FakeData(newDataPOT)));
179  else
180  ret.emplace("2019_beamPred", std::make_unique<Spectrum>(pred.Predict(calc).MockData(newDataPOT)));
181  auto & fakeData = ret.at("2019_beamPred");
182  fakeData->OverridePOT(FULL_RHC_EXPOSURE);
183  *sumPred += *fakeData;
184 
185  std::cout << " new " << (asimov ? "Asimov" : "fake") << " data = "
186  << fakeData->Integral(fakeData->POT(), nullptr, kPOT) << " events"
187  << " (" << (scaleFactor-1)/scaleFactor * FULL_RHC_EXPOSURE << " POT)" << std::endl;
188 
189  std::cout << " new total = " << sumPred->Integral(sumPred->POT(), nullptr, kPOT) << " events" << std::endl;
190 
191  return std::move(ret);
192  } // Get2019Prediction()
193 
194  } // namespace fakedata
195 } // namespace ana2019
196 
197 //--------------------------------------------------------------
198 
200 {
203 
204  std::cout << "\n\nThrowing mock data universes..." << std::endl;
205 
206  TFile outfile("fake_data_2019.root", "recreate");
207 
208  for (int univ = -1; univ < ana2019::fakedata::NUM_FAKEDATA_THROWS; univ++)
209  {
210  bool asimov = univ == -1;
211  std::cout << "\nUniv: " << (asimov ? "Asimov" : std::to_string(univ)) << std::endl;
212  std::cout << "============" << std::endl;
213  TDirectory * univDir = outfile.mkdir(asimov ? "asimov" : Form("univ_%d", univ));
214  for (auto analysis : {kNueAna2018, kNumuAna2018})
215  {
216  TDirectory * analysisDir = univDir->mkdir(analysis == kNueAna2018 ? "nue" : "numu");
217  for (auto beam : {kFHC, kRHC})
218  {
219  std::string beamStr = beam == kFHC ? "fhc" : "rhc";
220  TDirectory * beamDir = analysisDir->mkdir(beamStr.c_str());
221  if (analysis == kNueAna2018)
222  {
223  for (const auto & specPair : ana2019::fakedata::Get2019Prediction(&calc, analysis, beam, -1, asimov))
224  specPair.second->SaveTo(beamDir->mkdir(specPair.first.c_str()));
225  }
226  else if (analysis == kNumuAna2018)
227  {
228  for (auto q : {1, 2, 3, 4})
229  {
230  TDirectory * qDir = beamDir->mkdir(Form("q%d", q));
231  for (const auto & specPair : ana2019::fakedata::Get2019Prediction(&calc, analysis, beam, q, asimov))
232  specPair.second->SaveTo(qDir->mkdir(specPair.first.c_str()));
233  }
234  }
235  } // for (beam)
236  } // for (analysis)
237  } // for (univ)
238  outfile.Close();
239 
240  std::cout << "\n\ndone." << std::endl;
241 
242 }
const XML_Char * name
Definition: expat.h:151
Shared items for fake data studies for 2019 analysis.
const StorageType & Get(EAnaType2018 analysis, BeamType2018 beam, int quartileIfNumu=-1) const
Note: quartileIfNumu counts from 1, not 0.
osc::OscCalculatorDumb calc
Oscillation analysis framework, runs over CAF files outside of ART.
osc::IOscCalculatorAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
const TieredAnaObjStorage< Spectrum > kDataSpectra([](const std::string &beamStr){return GetNueData2018(beamStr);}, [](const std::string &beamStr){auto vec=GetNumuData2018(4, beamStr);std::vector< const Spectrum * > ret;ret.insert(ret.end(), vec.begin(), vec.end());return ret;})
General interface to oscillation calculators.
Definition: FwdDeclare.h:15
Optimized version of OscCalculatorPMNS.
Definition: StanTypedefs.h:28
void OverridePOT(double newpot)
DO NOT USE UNLESS YOU ARE 110% CERTAIN THERE ISN&#39;T A BETTER WAY!
Definition: Spectrum.h:269
void generate_fake_data()
const TieredAnaObjStorage< IPrediction > kPredictions([](const std::string &beamStr){return GetNuePrediction2018(beamStr=="fhc"?"combo":"prop", DefaultOscCalc(), false, beamStr, false);}, [](const std::string &beamStr){return GetNumuPredictions2018(4, false, beamStr);})
std::vector< Spectrum * > GetNumuData2018(const int nq=4, std::string beam="fhc")
void SetAna2018Calc(osc::IOscCalculatorAdjustable &calc)
Definition: fake_data.h:21
OStream cerr
Definition: OStream.cxx:7
Store the 2019 objects in tiered storage.
string filename
Definition: shutoffs.py:106
std::vector< const IPrediction * > GetNumuPredictions2018(const int nq=4, bool useSysts=true, std::string beam="fhc", bool GetFromUPS=false, ENu2018ExtrapType numuExtrap=kNuMu, bool minimizeMemory=false, bool NERSC=false)
virtual Spectrum Predict(osc::IOscCalculator *calc) const =0
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:37
const double kAna2018RHCPOT
Definition: Exposures.h:208
const double FULL_RHC_EXPOSURE
Definition: fake_data.h:17
Spectrum FakeData(double pot) const
Fake data is a MC spectrum scaled to the POT expected in the data.
Definition: Spectrum.cxx:819
std::map< std::string, std::unique_ptr< Spectrum > > Get2019Prediction(osc::IOscCalculator *calc, EAnaType2018 analysis, BeamType2018 beam, int quartileIfNumu=-1, bool asimov=false)
const IPrediction * GetNuePrediction2018(std::string decomp, osc::IOscCalculator *calc, bool corrSysts, std::string beam, bool isFakeData, bool GetFromUPS=false, bool minimizeMemory=false, bool NERSC=false)
Eigen::VectorXd vec
std::map< EAnaType2018, std::map< BeamType2018, std::vector< std::unique_ptr< const StorageType > > > > fStoredObjs
TieredAnaObjStorage(const std::function< const StorageType *(const std::string &beamStr)> &nueGetter, const std::function< std::vector< const StorageType * >(const std::string &beamStr)> &numuGetter)
OStream cout
Definition: OStream.cxx:6
const TieredAnaObjStorage< Spectrum > kCosmicSpectra([](const std::string &beamStr){std::string name;std::string anaDir="/nova/ana/nu_e_ana/Ana2018/";if(beamStr=="rhc") name="cosmic_spect_rhc";else name="cosmic_spect_fhc";auto spec=LoadFromFile< Spectrum >(anaDir+"/Predictions/cosmic/cosmic_prediction_real_data.root", name).release();return spec;}, [](const std::string &beamStr){std::string dir="/nova/ana/nu_mu_ana/Ana2018/Cosmics/";std::string filename=(beamStr=="rhc")?"cosmics_rhc__numu2018.root":"cosmics_fhc__numu2018.root";auto fcosm=TFile::Open((dir+filename).c_str());if(fcosm->IsZombie()){std::cerr<< "bad cosmics\n";exit(1);}std::vector< const Spectrum * > numu_cosmics;for(int i=1;i<=4;++i) numu_cosmics.emplace_back(new Spectrum(dynamic_cast< TH1 * >(fcosm->Get(Form("cosmics_q%d", i))), 0, (beamStr=="rhc")?kAna2018RHCLivetime:kAna2018FHCLivetime));return numu_cosmics;})
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
Spectrum MockData(double pot, int idx=0) const
Mock data is FakeData with Poisson fluctuations applied.
Definition: Spectrum.cxx:787
TDirectory * dir
Definition: macro.C:5
const int NUM_FAKEDATA_THROWS
Definition: fake_data.h:19
exit(0)
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
Spectrum * GetNueData2018(std::string beam)
const double kAna2018FHCLivetime
Definition: Exposures.h:213
FILE * outfile
Definition: dump_event.C:13
const double kAna2018RHCLivetime
Definition: Exposures.h:214