MakeSurface.C
Go to the documentation of this file.
1 /*
2  * MakeSurface.C
3  *
4  * Created on: May 7, 2018
5  *
6  * Author: J. Hewes <jhewes15@fnal.gov>
7  *
8  * Create contours for 2019 covmx analysis
9  */
10 
12 #include "CAFAna/Core/Progress.h"
13 #include "CAFAna/Core/Sample.h"
14 
16 
19 
20 #include "CAFAna/Analysis/CovMxSurface.h"
23 #include "CAFAna/Fit/Fit.h"
24 
26 
27 #include <ctime>
28 
29 using namespace ana;
30 
31 void MakeSurface(TString opt) {
32 
33  std::cout << "MakeSurface: Begin." << std::endl;
34 
35  DontAddDirectory guard;
36 
37  /////////////////////////////////////////////////////////////////////
38  // //
39  // ARGUMENT PARSING //
40  // //
41  /////////////////////////////////////////////////////////////////////
42 
43  std::vector<std::string> polarity = { "fhc", "rhc" };
44  std::cout << "MakeSurface: Getting Polarity Option." << std::endl;
45  std::vector<bool> use_polarity = GetOptionConfig(opt, polarity);
46  std::vector<std::string> detector = { "neardet", "fardet" };
47  std::cout << "MakeSurface: Getting Detector Option." << std::endl;
48  std::vector<bool> use_detector = GetOptionConfig(opt, detector);
49  std::cout << "MakeSurface: ...done" << std::endl;
50 
51  std::string optstring(opt.Data());
52 
53  // Get options for a particular systematic treatment.
54  // Expect:
55  // Single cm syst : _covmxsyst=name_
56  // Single pull syst : _pullsyst=name_
57  // All cm systs : _covmxsysts_
58  // All pull systs : _pullsyst=all_
59  // All cm, single pull swapped : _covmxsysts_pullsyst=name_
60 
61  // Single pull syst
62  std::string pull_syst = "none";
63  size_t pspos = optstring.find("pullsyst=");
64  if (pspos != std::string::npos) {
65  if (optstring.substr(pspos+8, 1) != "=")
66  assert(false and "The pullsyst option must take the form \"pullsyst=xxx\".");
67  size_t start = pspos+9;
68  size_t end = optstring.find("_", pspos);
69  if (end != std::string::npos)
70  pull_syst = optstring.substr(start, end-start);
71  else
72  pull_syst = std::stod(optstring.substr(start));
73  std::cout << "Treating syst " << pull_syst << " as a pull term." << std::endl;
74  }
75 
76  // All CM systs or All Pull systs
77  bool covmx_systs = CheckOption(opt, "covmxsysts");
78  if (pull_syst == "all" and covmx_systs) assert (false and
79  "Treating systematics as both pull terms and covariance matrix is not permitted.");
80 
81  // Single CM syst
82  std::string cm_syst = "none";
83  size_t cmspos = optstring.find("covmxsyst=");
84  if (cmspos != std::string::npos) {
85  if (optstring.substr(cmspos+9, 1) != "=")
86  assert(false and "The covmxsyst option must take the form \"covmxsyst=xxx\".");
87  if (covmx_systs)
88  assert(false and "All covmxsysts already, \"covmxsyst=xxx\" is redundant");
89  if (pull_syst!="none")
90  assert(false and "Are you sure you want to use a single CM syst with a single pull term?");
91  size_t start = cmspos+10;
92  size_t end = optstring.find("_", cmspos);
93  if (end != std::string::npos)
94  cm_syst = optstring.substr(start, end-start);
95  else
96  cm_syst = std::stod(optstring.substr(start));
97  std::cout << "Treating syst " << cm_syst << " as a cm term." << std::endl;
98  }
99 
100  // Load saved fake data, if requested
101  TDirectory* fakedata_dir = nullptr;
102  size_t upos = optstring.find("universe");
103  if (upos != std::string::npos) {
104  size_t universe;
105  size_t start = upos + 8;
106  size_t end = optstring.find("_", upos);
107  if (end != std::string::npos)
108  universe = std::stoi(optstring.substr(start, end - start));
109  else
110  universe = std::stoi(optstring.substr(start));
111 
112  std::cout << "Using saved fake data universe " << universe << std::endl;
113 
114  std::string fakedata_filename = FindCAFAnaDir() + "data/joint/fakedata.root";
115  TFile* fakedata_file = TFile::Open(fakedata_filename.c_str(), "read");
116  std::string dirname = "universe_" + std::to_string(universe);
117  fakedata_dir = (TDirectory*)fakedata_file->Get(dirname.c_str());
118  if (!fakedata_dir)
119  assert(false and "Unable to load fake data directory.");
120  }
121 
122  // Set delta m, if requested
123  double dm41 = -1;
124  size_t dmpos = optstring.find("dm41=");
125  if (dmpos != std::string::npos) {
126  if (optstring.substr(dmpos+4, 1) != "=")
127  throw std::runtime_error("The dm41 option must take the form \"dm41=xxx\".");
128  size_t start = dmpos+5;
129  size_t end = optstring.find("_", dmpos);
130  if (end != std::string::npos)
131  dm41 = std::stod(optstring.substr(start, end-start));
132  else
133  dm41 = std::stod(optstring.substr(start));
134  std::cout << "Using dm41 value " << dm41 << std::endl;
135  }
136 
137  // Don't profile delta24
138  bool prof_delta24 = not CheckOption(opt, "noprofdelta24");
139 
140  // Check what type of plot we're making
141  bool th24vsth34 = false;
142  bool th24vsdm41 = false;
143  bool th34vsdm41 = false;
144  if (CheckOption(opt, "th24vsth34")) th24vsth34 = true;
145  else if (CheckOption(opt, "th24vsdm41")) th24vsdm41 = true;
146  else if (CheckOption(opt, "th34vsdm41")) th34vsdm41 = true;
147  else throw std::runtime_error("The options string must contain plot type: \"th24vsth34\", \"th24vsdm41\" or \"th34vsdm41\".");
148 
149  // Make sure dm41 is set for th24 vs th34 case
150  if (th24vsth34 && dm41 == -1)
151  throw std::runtime_error("If you want to run a th24 vs th34 surface, you must specify a value for dm41.");
152 
153  /////////////////////////////////////////////////////////////////////
154  // //
155  // PREPARE INPUTS //
156  // //
157  /////////////////////////////////////////////////////////////////////
158 
159  // Oscillation calculator
161  SetAngles(calc);
162 
163  // Load predictions
164  std::vector<covmx::Sample> samples;
165  std::vector<IPrediction*> preds;
166  std::vector<const ISyst*> isysts;
167  for (size_t d = 0; d < 2; ++d) { // for detector
168  for (size_t p = 0; p < 2; ++p) { // for beam polarity
169  if (!use_polarity[p] || !use_detector[d]) continue;
170  samples.push_back(covmx::Sample(covmx::Analysis::kNC,
172  if(samples.size()==1)
173  isysts = GetSystsFromFile(samples.back());
174  std::cout << "MakeSurface: Loading Prediction..." << std::endl;
175  preds.push_back(LoadPrediction(samples.back(), pull_syst == "none" ? "nosysts" : "systs"));
176  std::cout << "MakeSurface: ...Loaded Prediction." << std::endl;
177  } // for beam polarity
178  } // for detector
179 
180  // Load cosmics
181  std::vector<Spectrum*> cosmic;
182  for (auto sample : samples){
183  std::cout << "MakeSurface: Loading Cosmic..." << std::endl;
184  cosmic.push_back(LoadCosmic(sample).release());
185  std::cout << "MakeSurface: ...Loaded Cosmic." << std::endl;
186  }
187 
188  // Load data
189  std::vector<Spectrum*> data;
190  if (fakedata_dir) {
191  for (auto sample : samples){
192  std::cout << "MakeSurface: Loading Fake Data..." << std::endl;
193  data.push_back(Spectrum::LoadFrom(fakedata_dir, sample.GetName().c_str())).release();
194  std::cout << "MakeSurface: ...Loaded Fake Data." << std::endl;
195  }
196  }
197  else {
198  for (size_t i = 0; i < preds.size(); ++i) {
199  data.push_back(new Spectrum(preds[i]->Predict(calc).FakeData(12e20)));
200  if(samples[i].detector == covmx::Detector::kFarDet)
201  data.back()->OverrideLivetime(kNus19SensLivetime);
202  if (cosmic[i]) *data[i] += *cosmic[i];
203  }
204  }
205 
206  // Load systematics
207  //std::vector<const ISyst*> isysts = GetSystsFromFile(samples[0]);
208 
209  // Pull term systematics
210  std::vector<const ISyst*> pull_systs = {};
211  SystConcat* sc = nullptr;
212  if (pull_syst == "all") {
213  // All pull systs, single detector
214  if (samples.size() == 1)
215  pull_systs = isysts;
216  // All pull systs, both detectors
217  else {
218  sc = new SystConcat(GetSystConcat(samples, isysts, preds));
219  pull_systs = sc->GetActiveSysts();
220  }
221  }
222  else if (pull_syst != "none") {
223  // Single pull syst, single detector
224  if (samples.size() == 1) {
225  for (const ISyst* syst : isysts) {
226  if (GetSystBaseName(syst->ShortName()) == pull_syst) {
227  pull_systs = { syst };
228  break;
229  }
230  }
231  if (pull_systs.empty()) {
232  std::ostringstream err;
233  err << "Couldn't find requested pull systematic " << pull_syst;
234  assert (false and err.str().c_str());
235  }
236  }
237  // Single pull syst, both detectors
238  else {
239  std::vector<SystConcat> all_scs = GetSystConcats(samples, preds, isysts);
240  for (auto it_sc : all_scs) {
241  if (GetSystBaseName(it_sc.GetActiveSysts()[0]->ShortName()) == pull_syst) {
242  sc = new SystConcat(it_sc);
243  pull_systs = sc->GetActiveSysts();
244  break;
245  }
246  }
247  }
248  }
249 
250  // Load covariance matrix
251  CovarianceMatrix* mx = nullptr;
252  if (covmx_systs || cm_syst!="none") {
253  std::vector<CovarianceMatrix*> CMs;
254  std::vector<Binning> bins;
255  for (auto d : data) bins.push_back(d->GetBinnings()[0]);
256  mx = new CovarianceMatrix(bins);
257  std::cout << std::endl;
258  for (CovarianceMatrix* it_mx : GetMatrices(samples)) {
259  std::string sysname = it_mx->GetName();
260  // if single CM syst configured, ignore all others
261  if (cm_syst!="none" && cm_syst!=sysname) continue;
262  // otherwise add the syst as long as it's not a pull term
263  if (GetSystBaseName(sysname) != pull_syst) {
264  std::cout << "Adding " << sysname
265  << " to covariance matrix." << std::endl;
266  mx->AddMatrix(it_mx);
267  }
268  else
269  std::cout << "Omitting pull term " << pull_syst
270  << " from matrix." << std::endl;
271  }
272  }
273 
274  /////////////////////////////////////////////////////////////////////
275  // //
276  // MAKE SURFACE //
277  // //
278  /////////////////////////////////////////////////////////////////////
279 
280  // Create experiment
281  OscCovMxExperiment* expt = nullptr;
282  if (covmx_systs||cm_syst!="none") expt = new OscCovMxExperiment(preds, mx, data, cosmic);
283  else expt = new OscCovMxExperiment(preds, data, cosmic);
284 
285  // Set Gaussian constraints, make MultiExperiment
286  std::vector<const IExperiment*> expts = GetConstraints();
287  expts.push_back(expt);
288  MultiExperiment multi_expt(expts);
289 
290  // Set fixed dm41 for th24 vs th34 case
291  if (th24vsth34){
292  calc->SetDm(4, dm41);
293  std::cout << "\n\n Set the dm41 \n\n" << std::endl;
294  }
295 
296  // Binning, fit vars and seeds
297  std::vector<const IFitVar*> fit_vars = { &kFitTheta23Sterile,
299  if (prof_delta24) fit_vars.push_back(&kFitDelta24InPiUnitsSterile);
300  else calc->SetDelta(2, 4, 0.5*M_PI);
301  std::map<const IFitVar*, std::vector<double> > seed_values;
302  seed_values[&kFitTheta23Sterile] = { kNus19Th23 };
303  seed_values[&kFitDmSq32Sterile] = { kNus19Dm32 };
304  if (prof_delta24) seed_values[&kFitDelta24InPiUnitsSterile] = { 0, 0.5, 1.5 };
305 
306  Binning* xbins = nullptr;
307  Binning* ybins = nullptr;
308 
309  std::vector<const IFitVar*> xy_vars;
310 
311  if (th24vsth34) {
312  xbins = new Binning(Binning::Simple(50, 0, 50));
313  ybins = new Binning(Binning::Simple(50, 0, 45));
314  xy_vars.push_back(&kFitTheta34InDegreesSterile);
315  xy_vars.push_back(&kFitTheta24InDegreesSterile);
316  }
317  else if (th24vsdm41 || th34vsdm41) {
318  xbins = new Binning(Binning::LogUniform(50, 1e-4, 1));
319  ybins = new Binning(Binning::LogUniform(50, 1e-2, 100));
320  if (th24vsdm41) {
321  xy_vars.push_back(&kFitSinSqTheta24Sterile);
322  xy_vars.push_back(&kFitDmSq41Sterile);
323  fit_vars.push_back(&kFitSinSqTheta34Sterile);
324  seed_values[&kFitSinSqTheta34Sterile] = { 0.01, 0.5, 0.99 };
325  }
326  if (th34vsdm41) {
327  xy_vars.push_back(&kFitSinSqTheta34Sterile);
328  xy_vars.push_back(&kFitDmSq41Sterile);
329  fit_vars.push_back(&kFitSinSqTheta24Sterile);
330  seed_values[&kFitSinSqTheta24Sterile] = { 0.01, 0.5, 0.99 };
331  }
332  }
333 
334  // Create surface
335  size_t part = 1275;
336  size_t njobs = 2500;
337  if (RunningOnGrid()){
338  part=JobNumber();
339  njobs=25;
340  }
341 
342  CovMxSurface surf(*xbins, *ybins, &multi_expt, calc, xy_vars[0], xy_vars[1],
343  fit_vars, seed_values, pull_systs, SystShifts::Nominal(), true, nullptr,
344  preds, data, cosmic, part, njobs);
345 
346  // Open output file and save surface
347  std::string outfile = optstring + ".root";
348  TFile* f = new TFile(outfile.c_str(), "recreate");
349  surf.SaveTo(f, optstring.c_str());
350  delete f;
351 
352 }
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 double kNus19Dm32
Definition: Utilities.h:28
const FitSinSqTheta24Sterile kFitSinSqTheta24Sterile
size_t JobNumber()
What&#39;s the process number for a grid job?
Definition: Utilities.cxx:438
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
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
void MakeSurface()
Definition: MakeSurface.C:31
double stod(const std::string &val)
bool RunningOnGrid()
Is this a grid (condor) job?
Definition: Utilities.cxx:359
const char * p
Definition: xmltok.h:285
void SetDelta(int i, int j, double delta)
Adapt the PMNS_Sterile calculator to standard interface.
const FitDmSq32Sterile kFitDmSq32Sterile
const FitTheta24InDegreesSterile kFitTheta24InDegreesSterile
const FitSinSqTheta34Sterile kFitSinSqTheta34Sterile
static SystShifts Nominal()
Definition: SystShifts.h:34
osc::OscCalcDumb calc
std::string FindCAFAnaDir()
Definition: Utilities.cxx:204
void SetAngles(osc::OscCalcSterile *calc)
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
#define M_PI
Definition: SbMath.h:34
void SaveTo(TDirectory *dir) const
Standard SaveTo implementation.
const XML_Char const XML_Char * data
Definition: expat.h:268
const FitDelta24InPiUnitsSterile kFitDelta24InPiUnitsSterile
TString part[npart]
Definition: Style.C:32
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.
std::vector< const IExperiment * > GetConstraints()
Gaussian constrains on atmospheric parameters.
Definition: Utilities.h:416
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
const int xbins
Definition: MakePlots.C:82
void FakeData()
Definition: rootlogon.C:156
std::vector< bool > GetOptionConfig(TString opt, std::vector< std::string > options)
Compare a single data spectrum to the MC + cosmics expectation.
Float_t d
Definition: plot.C:236
std::vector< CovarianceMatrix * > GetMatrices(std::vector< covmx::Sample > samples, bool cvmfs=true)
Definition: Utilities.h:350
std::vector< float > Spectrum
Definition: Constants.h:728
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 int ybins
const Binning bins
Definition: NumuCC_CPiBin.h:8
Combine multiple component experiments.
const double kNus19Th23
Definition: Utilities.h:33
void SetDm(int i, double dm)
IPrediction * LoadPrediction(std::string detector, bool rhc=false, std::string syst_type="all")
Function to load prediction object.
static Binning LogUniform(int n, double lo, double hi)
Definition: Binning.cxx:118
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)
const FitTheta34InDegreesSterile kFitTheta34InDegreesSterile
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
Float_t e
Definition: plot.C:35
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
const FitDmSq41Sterile kFitDmSq41Sterile
const double kNus19SensLivetime
Definition: Utilities.h:39
FILE * outfile
Definition: dump_event.C:13
std::vector< const ISyst * > GetSystsFromFile(covmx::Sample sample)
Get systematics for a given sample.
const FitTheta23Sterile kFitTheta23Sterile
surf
Definition: demo5.py:35
enum BeamMode string