BestFit.C
Go to the documentation of this file.
1 /*
2  * BestFit.C
3  *
4  * Created on: August 23, 2019
5  *
6  * Author: J. Hewes <jhewes15@fnal.gov>
7  *
8  * Find best fit for 2020 covmx analysis
9  */
10 
11 #include "CAFAna/Core/Progress.h"
12 #include "CAFAna/Core/Sample.h"
16 #include "CAFAna/Fit/Fit.h"
19 
21 
22 #include <ctime>
23 
24 using namespace ana;
25 
26 void BestFit(TString opt, std::string path = "/nova/data/users/jhewes15/nus20/fits") {
27 
28  try {
29 
30  DontAddDirectory guard;
31 
32  /////////////////////////////////////////////////////////////////////
33  // //
34  // ARGUMENT PARSING //
35  // //
36  /////////////////////////////////////////////////////////////////////
37 
38  std::string optstring(opt.Data());
39 
40  // Get pull and covariance matrix syst options
41  bool covmxSysts = CheckOption(opt, "covmxsysts");
42  std::string pullSyst = "none";
43  size_t pos = optstring.find("pullsyst=");
44  if (pos != std::string::npos) {
45  size_t start = pos + 9;
46  size_t end = optstring.find("_", pos);
47  pullSyst = optstring.substr(start, end-start);
48  std::cout << "Treating syst " << pullSyst << " as a pull term." << std::endl;
49  }
50  if (pullSyst == "all" and covmxSysts) assert (false and
51  "Treating systematics as both pull terms and covariance matrix is not permitted.");
52 
53  // Get samples, predictions and systematics
54  std::vector<covmx::Sample> samples = GetSamplesFromOptString(opt);
55  if (samples.size() == 0) {
56  std::cerr << "No valid samples found in option string! Exiting." << std::endl;
57  assert(false);
58  }
59  std::vector<std::vector<const ISyst*>> isysts;
60  for (covmx::Sample& sample : samples) {
61  sample.SetPrediction(LoadPrediction(sample, pullSyst == "none" ? "nosysts" : "systs", kUseCVMFS));
62  sample.SetPOT(sample.detector == covmx::kNearDet ? kAna2018SensitivityFHCNDPOT : kAna2018FHCPOT);
63  sample.SetLivetime(sample.detector == covmx::kNearDet ? -1 : kAna2018FHCLivetime);
64  isysts.push_back(LoadSysts(sample, kUseCVMFS));
65  }
66 
67  // Load saved fake data, if requested
68  std::string fakedata;
69  size_t universe = 999999;
70  pos = optstring.find("fakedata=");
71  if (pos != std::string::npos) {
72  size_t start = pos + 9;
73  if (optstring.substr(start, 6) == "asimov") {
74  fakedata = "asimov";
75  }
76  else if (optstring.substr(start, 5) == "stats") {
77  fakedata = "stats";
78  size_t end = optstring.find("_", pos);
79  universe = std::stoi(optstring.substr(start+5, end-(start+5)));
80  }
81  else if (optstring.substr(start, 5) == "systs") {
82  fakedata = "systs";
83  size_t end = optstring.find("_", pos);
84  universe = std::stoi(optstring.substr(start+5, end-(start+5)));
85  }
86  else {
87  assert(false and "Fake data type not recognised!");
88  }
89  }
90 
91  // Set delta m, if requested
92  double dm41 = -1;
93  pos = optstring.find("dm41=");
94  if (pos != std::string::npos) {
95  size_t start = pos+5;
96  size_t end = optstring.find("_", start);
97  dm41 = std::stod(optstring.substr(start, end-start));
98  std::cout << "Using dm41 value " << dm41 << std::endl;
99  }
100 
101  // Don't profile delta24
102  bool profDelta24 = not CheckOption(opt, "noprofdelta24");
103  bool profVars = not CheckOption(opt, "noprofvars");
104 
105  // Zoom y axis
106  std::tuple<bool, double, double> yrange = std::make_tuple(false, -1, -1);
107  pos = optstring.find("zoomyaxis");
108  if (pos != std::string::npos) {
109  // Using manually defined y axis range
110  std::get<0>(yrange) = true;
111  // Get lower bound
112  size_t start = pos+10;
113  size_t end = optstring.find("-", start);
114  // std::cout << "Lower bound: " << optstring.substr(start, end-start) << std::endl;
115  std::get<1>(yrange) = std::stod(optstring.substr(start, end-start));
116  // Get upper bound
117  start = end + 1;
118  end = optstring.find("_", start);
119  // std::cout << "Upper bound: " << optstring.substr(start, end-start) << std::endl;
120  std::get<2>(yrange) = std::stod(optstring.substr(start, end-start));
121  }
122 
123 
124  // Check what type of plot we're making
125  bool th24vsth34 = false;
126  bool th24vsdm41 = false;
127  bool th34vsdm41 = false;
128  if (CheckOption(opt, "th24vsth34")) th24vsth34 = true;
129  else if (CheckOption(opt, "th24vsdm41")) th24vsdm41 = true;
130  else if (CheckOption(opt, "th34vsdm41")) th34vsdm41 = true;
131  else throw std::runtime_error("The options string must contain plot type: \"th24vsth34\", \"th24vsdm41\" or \"th34vsdm41\".");
132 
133  // Make sure dm41 is set for th24 vs th34 case
134  if (th24vsth34 && dm41 == -1) {
135  throw std::runtime_error("If you want to run a th24 vs th34 surface, you must specify a value for dm41.");
136  }
137 
138  /////////////////////////////////////////////////////////////////////
139  // //
140  // PREPARE INPUTS //
141  // //
142  /////////////////////////////////////////////////////////////////////
143 
144  // Oscillation calculator
146  SetNus20Params(calc);
147  PrintOscParams(calc);
148 
149  // Load cosmics
150  std::vector<Spectrum*> cosmic;
151  for (auto sample : samples) {
152  cosmic.push_back(LoadCosmic(sample, kUseCVMFS).release());
153  }
154 
155  // Load data
156  std::vector<Spectrum*> data;
157  if (fakedata == "asimov") {
158  std::cout << "Using Asimov fake data." << std::endl;
159  for (size_t i = 0; i < samples.size(); ++i) {
160  data.push_back(new Spectrum(samples[i].GetPrediction()->Predict(calc).FakeData(samples[i].GetPOT())));
161  if(samples[i].detector == covmx::kFarDet)
162  data.back()->OverrideLivetime(samples[i].GetLivetime());
163  if (cosmic[i]) *data[i] += *cosmic[i];
164  }
165  }
166  else {
167  std::cout << "Loading universe " << universe << " of " << fakedata
168  << " fake data." << std::endl;
169  data = LoadFakeData(samples, universe, (fakedata=="systs"), kUseCVMFS);
170  }
171 
172  // Pull term systematics
173  std::vector<const ISyst*> pullSysts = {};
174  SystConcat* sc = nullptr;
175  if (pullSyst == "all") {
176  // All pull systs, single detector
177  if (samples.size() == 1)
178  pullSysts = isysts[0];
179  // All pull systs, both detectors
180  else {
181  sc = new SystConcat(GetSystConcat(samples, isysts));
182  pullSysts = sc->GetActiveSysts();
183  }
184  }
185  else if (pullSyst != "none") {
186  // Single pull syst, single detector
187  if (samples.size() == 1) {
188  for (std::vector<const ISyst*> detSysts : isysts) {
189  for (const ISyst* syst : detSysts) {
190  if (GetSystBaseName(syst->ShortName()) == pullSyst) {
191  pullSysts = { syst };
192  break;
193  }
194  }
195  }
196  if (pullSysts.empty()) {
197  std::ostringstream err;
198  err << "Couldn't find requested pull systematic " << pullSyst;
199  assert (false and err.str().c_str());
200  }
201  }
202  // Single pull syst, both detectors
203  else {
204  std::vector<SystConcat> allSCs = GetSystConcats(samples, isysts);
205  for (auto it_sc : allSCs) {
206  if (GetSystBaseName(it_sc.GetActiveSysts()[0]->ShortName()) == pullSyst) {
207  sc = new SystConcat(it_sc);
208  pullSysts = sc->GetActiveSysts();
209  break;
210  }
211  }
212  }
213  }
214 
215  // Load covariance matrix
216  covmx::CovarianceMatrix* mx = nullptr;
217  if (covmxSysts) {
218  for (size_t i = 0; i < data.size(); ++i) {
219  samples[i].SetBinning(data[i]->GetBinnings()[0]);
220  }
221  mx = new covmx::CovarianceMatrix(samples);
222  for (covmx::CovarianceMatrix* itMx : GetMatrices(samples, kUseCVMFS)) {
223  if (GetSystBaseName(itMx->GetName()) != pullSyst) {
224  mx->AddMatrix(itMx);
225  }
226  else
227  std::cout << "Omitting pull term " << pullSyst
228  << " from matrix." << std::endl;
229  }
230  }
231 
232  /////////////////////////////////////////////////////////////////////
233  // //
234  // MAKE SURFACE //
235  // //
236  /////////////////////////////////////////////////////////////////////
237 
238  // Create experiment
239  OscCovMxExperiment* expt = nullptr;
240  if (covmxSysts) expt = new OscCovMxExperiment(samples, mx, data, cosmic);
241  else expt = new OscCovMxExperiment(samples, data, cosmic);
242  if (samples.size() > 1 and sc) expt->ConfigureSysts(sc);
243 
244  // Set Gaussian constraints, make MultiExperiment
245  std::vector<const IChiSqExperiment*> expts = GetConstraints();
246  expts.push_back(expt);
247  MultiExperiment multiExpt(expts);
248 
249  // Binning, fit vars and seeds
250  std::string inFileName = path + "/" + optstring + ".root";
251  TFile* inFile = TFile::Open(inFileName.c_str(), "read");
252  CovMxSurface* surf = CovMxSurface::LoadFrom(inFile).release();
253  delete inFile;
254  std::vector<const IFitVar*> fitVars = { &kFitTheta23Sterile,
256  if (profDelta24) fitVars.push_back(&kFitDelta24InPiUnitsSterile);
257  else calc->SetDelta(2, 4, 0.5);
258 
259  const IFitVar* xvar = nullptr;
260  const IFitVar* yvar = nullptr;
261 
262  if (th24vsth34) {
265  }
266  else if (th24vsdm41) {
267  xvar = &kFitSinSqTheta24Sterile;
268  yvar = &kFitDmSq41Sterile;
269  fitVars.push_back(&kFitSinSqTheta34Sterile);
270  }
271  else if (th34vsdm41) {
272  xvar = &kFitSinSqTheta34Sterile;
273  yvar = &kFitDmSq41Sterile;
274  fitVars.push_back(&kFitSinSqTheta24Sterile);
275  }
276 
277  // Override fit variables if requested
278  if (not profVars) {
279  fitVars.clear();
280  }
281 
282  surf->FindBestFit(&multiExpt, calc, xvar, yvar, fitVars, pullSysts);
283 
284  // Open output file and save surface
285  std::string outfile = optstring + "_bestfit.root";
286  TFile* f = new TFile(outfile.c_str(), "recreate");
287  surf->SaveTo(f->mkdir(optstring.c_str()));
288  delete f;
289 
290  } catch(const std::exception& e) {
291  std::cout << std::endl << "Exception occurred! " << e.what() << std::endl;
292  }
293 
294 }
std::vector< SystConcat > GetSystConcats(std::vector< covmx::Sample > samples, std::vector< IPrediction * > preds, std::vector< const ISyst * > systs)
Get correctly concatenated systs for ND and FD.
Definition: Utilities.h:294
const FitSinSqTheta24Sterile kFitSinSqTheta24Sterile
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::string GetSystBaseName(std::string name)
Get the last part of a systematic&#39;s ShortName.
Definition: Utilities.h:232
double stod(const std::string &val)
void SetDelta(int i, int j, double delta)
Adapt the PMNS_Sterile calculator to standard interface.
const FitDmSq32Sterile kFitDmSq32Sterile
OStream cerr
Definition: OStream.cxx:7
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
void AddMatrix(CovarianceMatrix *gen)
const double kAna2018SensitivityFHCNDPOT
Definition: Exposures.h:210
void PrintOscParams(osc::OscCalcSterile *calc)
const FitTheta24InDegreesSterile kFitTheta24InDegreesSterile
const FitSinSqTheta34Sterile kFitSinSqTheta34Sterile
osc::OscCalcDumb calc
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
void BestFit(TString opt, std::string path="/nova/data/users/jhewes15/nus20/fits")
Definition: BestFit.C:26
const XML_Char const XML_Char * data
Definition: expat.h:268
const FitDelta24InPiUnitsSterile kFitDelta24InPiUnitsSterile
ifstream inFile
Definition: AnaPlotMaker.h:34
osc::OscCalcSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
Definition: Calcs.cxx:97
expt
Definition: demo5.py:34
Sum up livetimes from individual cosmic triggers.
void SetNus20Params(osc::OscCalcSterile *calc)
Definition: Utilities.h:682
std::vector< const IExperiment * > GetConstraints()
Gaussian constrains on atmospheric parameters.
Definition: Utilities.h:416
Compare a single data spectrum to the MC + cosmics expectation.
inFileName
if we get here, we&#39;re doing the base definitions:
Definition: mkDefs.py:168
std::vector< CovarianceMatrix * > GetMatrices(std::vector< covmx::Sample > samples, bool cvmfs=true)
Definition: Utilities.h:350
std::vector< float > Spectrum
Definition: Constants.h:610
SystConcat GetSystConcat(std::vector< covmx::Sample > samples, std::vector< const ISyst * > systs, std::vector< IPrediction * > preds)
Get correctly concatenated systs for ND and FD.
Definition: Utilities.h:249
bool CheckOption(TString opts, TString opt)
OStream cout
Definition: OStream.cxx:6
const std::string path
Definition: plot_BEN.C:43
Combine multiple component experiments.
IPrediction * LoadPrediction(std::string detector, bool rhc=false, std::string syst_type="all")
Function to load prediction object.
std::vector< covmx::Sample > GetSamplesFromOptString(TString optString)
Function to take an option TString and return a vector of associated covmx::Samples.
Definition: Utilities.h:384
void SetBinning(std::vector< Sample > samples)
std::unique_ptr< Spectrum > LoadCosmic(covmx::Sample sample, bool cvmfs=true)
Get cosmics for a given sample.
Definition: Utilities.h:145
assert(nhit_max >=nhit_nbins)
std::vector< const ISyst * > LoadSysts(covmx::Sample sample)
Get systematics for a given sample.
Definition: Utilities.h:524
const FitTheta34InDegreesSterile kFitTheta34InDegreesSterile
Interface definition for fittable variables.
Definition: IFitVar.h:16
double GetPOT(bool isFHC=true)
void FakeData()
Definition: rootlogon.C:156
std::map< std::string, ana::Spectrum > LoadFakeData(const std::string &fakeDataFilename, const std::string &path)
static std::unique_ptr< CovMxSurface > LoadFrom(TDirectory *dir)
Standard LoadFrom implementation.
IPrediction * GetPrediction(SystPredType predType, const HistAxis *axis, const HistAxis *numuAxis, const Cut *FDCut, const Cut *NDCut, const Cut *numuCut, const SystShifts *shiftData, const SystShifts *shiftMC, const Var *weight, Loaders *loaders, Loaders *loadersND)
Definition: SystMaker.cxx:44
const double kAna2018FHCPOT
Definition: Exposures.h:207
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
Float_t e
Definition: plot.C:35
const FitDmSq41Sterile kFitDmSq41Sterile
const double kAna2018FHCLivetime
Definition: Exposures.h:213
FILE * outfile
Definition: dump_event.C:13
A class for generating a covariance matrices as a function of oscillation parameters and systematics ...
const FitTheta23Sterile kFitTheta23Sterile
surf
Definition: demo5.py:35
enum BeamMode string