make_fc_nus_surfs_nersc_2019.C
Go to the documentation of this file.
1 #include "TFile.h"
2 #include "TH2.h"
3 #include "TRandom3.h"
4 #include "TDirectory.h"
5 #include "TH1.h"
6 
7 #include <ctime>
8 #include <fstream>
9 #include <sstream>
10 
12 
13 
14 
16 #include "CAFAna/Core/Progress.h"
17 #include "CAFAna/Core/Sample.h"
18 
20 
24 
25 #include "CAFAna/Analysis/CovMxSurface.h"
28 #include "CAFAna/Fit/Fit.h"
29 
30 #include "CAFAna/FC/FCCollection.h"
31 #include "CAFAna/FC/FCPoint.h"
32 #include "CAFAna/FC/FCSurface.h"
33 
35 
36 using namespace ana;
37 
38 std::vector<double> LogBins(int nbins,double xlo, double xhi)
39 {
40  double log_spacing = exp( (log(xhi) - log(xlo))/(nbins) );
41  std::vector<double> binning;
42  for (int i = 0; i <= nbins; ++i) {
43  if (i == 0) binning.push_back(xlo);
44  else binning.push_back(log_spacing*binning[i-1]);
45  }
46  return binning;
47 }
48 
49 std::vector<double> LinBins(int nbins,double xlo, double xhi)
50 {
51  double lin_spacing = (xhi - xlo)/(nbins);
52  std::vector<double> binning;
53  for (int i = 0; i <= nbins; ++i) {
54  if (i == 0) binning.push_back(xlo);
55  else binning.push_back(lin_spacing+binning[i-1]);
56  }
57  return binning;
58 }
59 
60 void make_fc_nus_surfs_nersc_2019(int NPts = 1, int startidx = 0, int endidx = 0, bool nh = true,
61  int N = 1, int jid = 1, int npts_flush = 1,
62  std::string plot = "th24vsth34_UO" ,
63  TString beam = "neardetfardetFHCcovmx", bool corrSysts = true,
64  bool readFromFile = true, double max_time_in_sec = 3600)
65 {
66  //clock starts
67  std::clock_t tstart = std::clock();
68 
69  DontAddDirectory guard;
70 
71  FCCollection fccol;
72 
73  std::vector<std::string> polarity = { "fhc", "rhc" };
74  std::vector<bool> use_polarity = GetOptionConfig(beam, polarity);
75  std::vector<std::string> detector = { "neardet", "fardet" };
76  std::vector<bool> use_detector = GetOptionConfig(beam, detector);
77 
78  bool covmx = beam.Contains("covmx", TString::kIgnoreCase);
79  bool extrap = beam.Contains("extrap", TString::kIgnoreCase);
80  assert(covmx xor extrap);
81 
82  if (covmx && use_polarity[1])
83  throw std::runtime_error("RHC covariance matrix not supported yet.");
84 
85  if (extrap && use_polarity[0])
86  throw std::runtime_error("Extrapolated fit not supported for FHC.");
87 
88  std::cout << "beam:: " << beam << std::endl;
89  std::cout << "systs:: " << corrSysts << std::endl;
90  if (covmx && corrSysts) {
91  corrSysts = false;
92  std::cout << "Using covariance matrix, so disabling traditional systematics." << std::endl;
93  }
94 
95  bool upperOctant = false;
96  if(plot == "th24vsth34_UO" ||
97  plot == "th24_UO" ||
98  plot == "th34_UO" ||
99  plot == "dm41vsth24_UO" ||
100  plot == "dm41vsth34_UO")
101  {
102  upperOctant = true;
103  }
104 
105  bool th24 = false;
106  bool th34 = false;
107  bool th24vsth34 = false;
108  bool dm41vsth24 = false;
109  bool dm41vsth34 = false;
110  if(plot == "th24_UO" || plot == "th24_LO"){
111  th24 = true;
112  }
113  else if(plot == "th34_UO" || plot == "th34_LO"){
114  th34 = true;
115  }
116  else if(plot == "th24vsth34_UO" || plot == "th24vsth34_LO") th24vsth34 = true;
117  else if(plot == "dm41vsth24_UO" || plot == "dm41vsth24_LO") dm41vsth24 = true;
118  else if(plot == "dm41vsth34_UO" || plot == "dm41vsth34_LO") dm41vsth34 = true;
119 
120  beam += "extrap";
121  if(corrSysts) beam += "syst";
122 
123  // Get POT & livetimes
124  //double nd_pot = use_polarity[1] ? kAna2018SensitivityRHCNDPOT : kAna2018SensitivityFHCNDPOT;
125  double fd_pot = use_polarity[1] ? kAna2018RHCPOT : kAna2018FHCPOT;
126  double fd_livetime = use_polarity[1] ? kAna2018RHCLivetime : kAna2018FHCLivetime;
127 
128  IPrediction* pred = nullptr;
129  Spectrum* cosmic_ext = nullptr;
130 
131  std::vector<covmx::Sample> samples;
132  std::vector<IPrediction*> preds;
133  std::vector<Spectrum*> cosmic;
134  std::vector<Spectrum*> data;
135  std::vector<Spectrum*>::iterator iter;
136 
137  CovarianceMatrix* mx = nullptr;
138 
139  std::string pull_syst = "none";
140  std::vector<const ISyst*> systs = {};
141  SystConcat* sc = nullptr;
142 
143  // OscCalc for Sterile 3+1 Model
145  SetAngles(calc);
146 
147  if (covmx) {
148 
149  for (size_t d = 0; d < 2; ++d) { // for detector
150  for (size_t p = 0; p < 2; ++p) { // for beam polarity
151  if (!use_polarity[p] || !use_detector[d]) continue;
152  samples.push_back(covmx::Sample(covmx::Analysis::kNC,
154  preds.push_back(LoadPrediction(samples.back(), pull_syst == "none" ? "nosysts" : "systs"));
155  } // for beam polarity
156  } // for detector
157 
158  for (auto sample : samples)
159  cosmic.push_back(LoadCosmic(sample).release());
160 
161  std::vector<const ISyst*> isysts = GetSystsFromFile(samples[0]);
162 
163  if (pull_syst == "all") {
164  // All pull systs, single detector
165  if (samples.size() == 1)
166  systs = isysts;
167  // All pull systs, both detectors
168  else {
169  sc = new SystConcat(GetSystConcat(samples, isysts, preds));
170  systs = sc->GetActiveSysts();
171  }
172  }
173  else if (pull_syst != "none") {
174  // Single pull syst, single detector
175  if (samples.size() == 1) {
176  for (const ISyst* syst : isysts) {
177  if (GetSystBaseName(syst->ShortName()) == pull_syst) {
178  systs = { syst };
179  break;
180  }
181  }
182  if (systs.empty()) {
183  std::ostringstream err;
184  err << "Couldn't find requested pull systematic " << pull_syst;
185  assert (false and err.str().c_str());
186  }
187  }
188  // Single pull syst, both detectors
189  else {
190  std::vector<SystConcat> all_scs = GetSystConcats(samples, preds, isysts);
191  for (auto it_sc : all_scs) {
192  if (GetSystBaseName(it_sc.GetActiveSysts()[0]->ShortName()) == pull_syst) {
193  sc = new SystConcat(it_sc);
194  systs = sc->GetActiveSysts();
195  break;
196  }
197  }
198  }
199  }
200 
201  for (size_t i = 0; i < preds.size(); ++i) {
202  data.push_back(new Spectrum(preds[i]->Predict(calc).FakeData(12e20)));
203  if(samples[i].detector == covmx::Detector::kFarDet)
204  data.back()->OverrideLivetime(kNus19SensLivetime);
205  if (cosmic[i]) *data[i] += *cosmic[i];
206  }
207 
208  std::vector<Binning> bins;
209  for (auto d : data) bins.push_back(d->GetBinnings()[0]);
210  mx = new CovarianceMatrix(bins);
211  std::cout << std::endl;
212  for (CovarianceMatrix* it_mx : GetMatrices(samples)) {
213  if (GetSystBaseName(it_mx->GetName()) != pull_syst) {
214  std::cout << "Adding " << it_mx->GetName()
215  << " to covariance matrix." << std::endl;
216  mx->AddMatrix(it_mx);
217  }
218  else
219  std::cout << "Omitting pull term " << pull_syst
220  << " from matrix." << std::endl;
221  }
222 
223  } // if covmx
224  else {
225 
226  // load prediction for Extrap fit
227  pred = LoadExtrapPrediction(beam);
228 
229  // load cosmics for Extrap fit
230  //cosmic = new Spectrum(LoadCosmicSpectrum(rhc));
231  cosmic_ext = new Spectrum(LoadCosmicSpectrum(polarity[1]));
232  cosmic_ext->Scale(fd_livetime/cosmic_ext->Livetime());
233  cosmic_ext->OverridePOT(fd_pot);
234 
235  if (corrSysts)
236  //systs = rhc ? GetNus18RHCExtrapSysts() : GetNus18FHCFDSysts("all");
237  systs = LoadISysts();
238 
239  }
240 
241  std::cout << "Number of Systs:: " << systs.size() << std::endl;
242 
243  // Binnings
244  std::vector<double> th24Bins;
245  std::vector<double> th34Bins;
246  std::vector<double> dm2Bins;
247 
248  if(th24){
249  th24Bins = LinBins(180, 0., 45.);
250  }
251  else if(th34) {
252  th34Bins = LinBins(360, 0., 90.);
253  }
254  else if (th24vsth34) {
255  th24Bins = LinBins(45, 0., 45.);
256  th34Bins = LinBins(90, 0., 90.);
257  }
258  else if (dm41vsth24) {
259  th24Bins = LogBins(50, 1.0e-6, 1);
260  dm2Bins = LogBins(50, 1.0e-2, 100);
261  }
262  else if (dm41vsth34) {
263  th34Bins = LogBins(50, 1.0e-6, 1);
264  dm2Bins = LogBins(50, 1.0e-2, 100);
265  }
266 
267  double binxlowedge = 0.;
268  double binxuppedge = 0.;
269  double binylowedge = 0.;
270  double binyuppedge = 0.;
271  double binxcenter = 0.;
272  double binycenter = 0.;
273 
274  int TotBins = 4050;
275  int bins[TotBins];
276 
277  std::ifstream binlist;
278  std::string path = "./";
279  std::string binlist_fname = path+"binlist.txt";
280  binlist.open(binlist_fname);
281  int j = 0;
282  // check if the binlist exists
283  if(!binlist.good()){
284  std::cout << "Could not find "<< binlist_fname << std::endl;
285  exit(1);
286  }
287  if(readFromFile){
288  do
289  {
290  binlist >> bins[j];
291  j++;
292  }while(!binlist.eof());
293  binlist.close();
294  }
295  else for (int i = 0; i < TotBins; i++) bins[i]=i;
296 
297  int ibin = bins[startidx];
298  int fbin = bins[endidx];
299 
300  std::string binname = "_bins_"+std::to_string(ibin) + "_to_" + std::to_string(fbin);
301  std::string hierStr = nh ? "NH" : "IH";
302  std::string number = std::to_string(N);
303  std::string fcout = "FCCol_nus18_"+plot+"_"+hierStr+"_";
304  fcout += use_polarity[0] ? "FHC":"RHC";
305  fcout += corrSysts ? "_systs":"_stats";
306  fcout += binname+"_v"+number+'_'+std::to_string(jid)+".root";
307 
308  // finished loading
309  std::clock_t tload = std::clock();
310  double loadtime = ( tload - tstart ) / (double) CLOCKS_PER_SEC;
311 
312  TRandom3 rnd(0); // zero seeds randomly
313 
314  for (int i = startidx; i <= endidx; i++) {
315 
316  int iter_bin = bins[i];
317 
318  if(th24){
319  binxlowedge = th34Bins.front();
320  binxuppedge = th34Bins.back();
321  binylowedge = th24Bins[iter_bin];
322  binyuppedge = th24Bins[iter_bin+1];
323  binxcenter = 0.5*(binxlowedge+binxuppedge);
324  binycenter = 0.5*(binylowedge+binyuppedge);
325  }
326  else if(th34){
327  binxlowedge = th34Bins[iter_bin];
328  binxuppedge = th34Bins[iter_bin+1];
329  binylowedge = th24Bins.front();
330  binyuppedge = th24Bins.back();
331  binxcenter = 0.5*(binxlowedge+binxuppedge);
332  binycenter = 0.5*(binylowedge+binyuppedge);
333  }
334  else if (th24vsth34) {
335  int x_bins = th34Bins.size() - 1;
336  int x_idx = iter_bin % x_bins;
337  int y_idx = iter_bin / x_bins;
338  binxlowedge = th34Bins[x_idx];
339  binxuppedge = th34Bins[x_idx+1];
340  binylowedge = th24Bins[y_idx];
341  binyuppedge = th24Bins[y_idx+1];
342  binxcenter = 0.5*(binxlowedge+binxuppedge);
343  binycenter = 0.5*(binylowedge+binyuppedge);
344  }
345  else if (dm41vsth24) {
346  int x_bins = th24Bins.size() - 1;
347  int x_idx = iter_bin % x_bins;
348  int y_idx = iter_bin / x_bins;
349  binxlowedge = th24Bins[x_idx];
350  binxuppedge = th24Bins[x_idx+1];
351  binylowedge = dm2Bins[y_idx];
352  binyuppedge = dm2Bins[y_idx+1];
353  binxcenter = 0.5*(binxlowedge+binxuppedge);
354  binycenter = 0.5*(binylowedge+binyuppedge);
355  }
356  else if (dm41vsth34) {
357  int x_bins = th34Bins.size() - 1;
358  int x_idx = iter_bin % x_bins;
359  int y_idx = iter_bin / x_bins;
360  binxlowedge = th34Bins[x_idx];
361  binxuppedge = th34Bins[x_idx+1];
362  binylowedge = dm2Bins[y_idx];
363  binyuppedge = dm2Bins[y_idx+1];
364  binxcenter = 0.5*(binxlowedge+binxuppedge);
365  binycenter = 0.5*(binylowedge+binyuppedge);
366  }
367 
368  std::cout << "Bin X (low-center-up): " << binxlowedge << "-" << binxcenter << "-" << binxuppedge << std::endl;
369  std::cout << "Bin Y (low-center-up): " << binylowedge << "-" << binycenter << "-" << binyuppedge << std::endl;
370 
371  double profDelta13_PiUnits = 1.46;
372  double profssTh13 = pow( sin( kNus19Th13 ) , 2);
373  double profssTh23 = 0.470;
374  double profDmSq = kNus19Dm21 + kNus19Dm32;
375  double profDelta24_PiUnits = 1.;
376  if(upperOctant){
377  profDelta13_PiUnits = kNus19Delta13/M_PI;
378  profssTh23 = pow( sin( kNus19Th23 ) , 2);
379  }
380  double profTh24 = 0.;
381  double profTh34 = 0.;
382  double profDm41 = 0.;
383  if(th24 || th34 || th24vsth34) {
384  profTh24 = binycenter;
385  profTh34 = binxcenter;
386  profDm41 = 0.5;
387  }
388  else if(dm41vsth24) {
389  profTh24 = binxcenter;
390  profTh34 = 0.5;
391  profDm41 = binycenter;
392  }
393  else if(dm41vsth34) {
394  profTh24 = 0.5;
395  profTh34 = binxcenter;
396  profDm41 = binycenter;
397  }
398 
399  // Vector of parameters to fit: #theta_{23}, #theta_{34}, #theta_{24}, #deltam_{32}
400  std::vector<const IFitVar*> prof_vars;
401  prof_vars.push_back(&kFitSinSqTheta23Sterile); // to be constrained to PDG
402  prof_vars.push_back(&kFitDmSq32Sterile); // to be constrained to PDG
403  prof_vars.push_back(&kFitDelta24InPiUnitsSterile);
404 
405  if(th24) {
406  prof_vars.push_back(&kFitTheta34InDegreesSterile);
407  }
408  else if(th34) {
409  prof_vars.push_back(&kFitTheta24InDegreesSterile);
410  }
411  else if(dm41vsth24) {
412  prof_vars.push_back(&kFitSinSqTheta34Sterile);
413  }
414  else if(dm41vsth34) {
415  prof_vars.push_back(&kFitSinSqTheta24Sterile);
416  }
417 
418  std::vector<const IFitVar*> fitvars = prof_vars;
419 
420  if(th24) {
421  fitvars.push_back(&kFitTheta24InDegreesSterile);
422  }
423  else if(th34) {
424  fitvars.push_back(&kFitTheta34InDegreesSterile);
425  }
426  else if(th24vsth34) {
427  fitvars.push_back(&kFitTheta24InDegreesSterile);
428  fitvars.push_back(&kFitTheta34InDegreesSterile);
429  }
430  else if(dm41vsth24) {
431  fitvars.push_back(&kFitSinSqTheta24Sterile);
432  fitvars.push_back(&kFitDmSq41Sterile);
433  }
434  else if(dm41vsth34) {
435  fitvars.push_back(&kFitSinSqTheta34Sterile);
436  fitvars.push_back(&kFitDmSq41Sterile);
437  }
438 
439  for(int iTrial = 0; iTrial < NPts; ++iTrial){
440 
441  std::clock_t tistart = std::clock();
442  if (max_time_in_sec > 0.0){
443  double elapstime = ( tistart - tstart ) / (double) CLOCKS_PER_SEC;
444  if (elapstime >= max_time_in_sec) break;
445  }
446 
447  SetAngles(calc);
448 
449  std::cout << "On trial: " << iTrial << std::endl;
450 
451  auto trueX = rnd.Uniform(binxlowedge, binxuppedge);
452  auto trueY = rnd.Uniform(binylowedge, binyuppedge);
453  auto trueDelta24 = rnd.Uniform(0.,2.);
454 
455  double trueTh24 = 0.;
456  double trueTh34 = 0.;
457  double trueDm41 = 0.;
458 
459  int mock_seed = 0;
460  std::cout << "Mock Seed: " << mock_seed << std::endl;
461 
462  if(th24 || th34 || th24vsth34) {
463  trueTh24 = trueY;
464  trueTh34 = trueX;
465  trueDm41 = 0.5;
466  kFitTheta34InDegreesSterile.SetValue(calc, trueTh34);
467  kFitTheta24InDegreesSterile.SetValue(calc, trueTh24);
468  }
469  else if(dm41vsth24) {
470  trueDm41 = trueY;
471  trueTh24 = trueX;
472  trueTh34 = rnd.Uniform(0.,1.);
473  kFitSinSqTheta24Sterile.SetValue(calc, trueTh24);
474  kFitSinSqTheta34Sterile.SetValue(calc, trueTh34);
475  }
476  else if(dm41vsth34) {
477  trueDm41 = trueY;
478  trueTh34 = trueX;
479  trueTh24 = rnd.Uniform(0.,1.);
480  kFitSinSqTheta24Sterile.SetValue(calc, trueTh24);
481  kFitSinSqTheta34Sterile.SetValue(calc, trueTh34);
482  }
483  kFitDelta24InPiUnitsSterile.SetValue(calc, trueDelta24);
484  kFitDmSq41Sterile.SetValue(calc, trueDm41);
485 
486  // Generate Mock Data
487  Spectrum* mock_spec = nullptr;
488 
489  for (iter=data.begin(); iter != data.end(); ++iter){
490  delete *iter;
491  }
492  data.clear();
493 
494  // For CovMx Fit
495  if (covmx) {
496 
497  for (size_t i = 0; i < preds.size(); ++i) {
498  data.push_back(new Spectrum(preds[i]->Predict(calc).MockData(12e20)));
499  if(samples[i].detector == covmx::Detector::kFarDet)
500  data.back()->OverrideLivetime(kNus19SensLivetime);
501  if (cosmic[i]) *data[i] += *cosmic[i];
502  }
503 
504  } // if covmx
505  // For Extrap Fit
506  else {
507 
508  Spectrum prediction = pred->Predict(calc);
509  cosmic_ext->Scale(fd_livetime/cosmic_ext->Livetime());
510  cosmic_ext->OverridePOT(fd_pot);
511  Spectrum cos_spec = cosmic_ext->MockData(fd_pot);
512  mock_spec = new Spectrum(prediction.MockData(fd_pot) + cos_spec);
513 
514  } // if extrap
515 
516  const IExperiment* expt;
517  if (covmx) {
518  if(corrSysts) expt = new OscCovMxExperiment(preds, mx, data, cosmic);
519  else expt = new OscCovMxExperiment(preds, data, cosmic);
520  }
521  else {
522  expt = new SingleSampleExperiment(pred, *mock_spec, *cosmic_ext);
523  }
524 
525  std::vector<const IExperiment*> expts = GetConstraints();
526  expts.push_back(expt);
527  MultiExperiment multi_expt(expts);
528 
529  // -------- TRUE FIT -------- //
530  MinuitFitter fittrue(&multi_expt, prof_vars, systs);
531 
532  //double trueChi = fittrue.Fit(calc)->EvalMetricVal();
533  double trueChi = fittrue.Fit(calc, IFitter::kQuiet)->EvalMetricVal();
534  double trueDelta13 = kFitDelta13InPiUnitsSterile.GetValue(calc);
535  double truessTh13 = kFitSinSqTheta13Sterile.GetValue(calc);
536  double truessTh23 = kFitSinSqTheta23Sterile.GetValue(calc);
537  double trueDmSq = kFitDmSq32Sterile.GetValue(calc);
538  trueDelta24 = kFitDelta24InPiUnitsSterile.GetValue(calc);
539 
540  if(th24) {
541  trueTh34 = kFitTheta34InDegreesSterile.GetValue(calc);
542  }
543  else if(th34) {
544  trueTh24 = kFitTheta24InDegreesSterile.GetValue(calc);
545  }
546  else if(dm41vsth24) {
547  trueTh34 = kFitSinSqTheta34Sterile.GetValue(calc);
548  }
549  else if(dm41vsth34) {
550  trueTh24 = kFitSinSqTheta24Sterile.GetValue(calc);
551  }
552 
553  // -------- BEST FIT -------- //
554  MinuitFitter fit(&multi_expt, fitvars, systs);
555 
556  double bestChi = 1E20;
557  double bestDelta13 = 0;
558  double bestssTh13 = 0;
559  double bestssTh23 = 0;
560  double bestDmSq = 0;
561  double bestDelta24 = 0;
562  double bestTh24 = 0;
563  double bestTh34 = 0;
564  double bestDm41 = 0;
565 
566  if(th24 || th34 || th24vsth34) {
567  for (double th34_seed: {15, 30, 45, 60, 75}) {
568  for (double th24_seed: {9, 18, 27, 36}) {
569  std::vector<double> th34_seed_vec; th34_seed_vec.push_back(th34_seed);
570  std::vector<double> th24_seed_vec; th24_seed_vec.push_back(th24_seed);
571  std::map<const IFitVar*, std::vector<double> > seedPts;
572  seedPts[&kFitTheta34InDegreesSterile] = th34_seed_vec;
573  seedPts[&kFitTheta24InDegreesSterile] = th24_seed_vec;
574 
575  SystShifts auxShifts = SystShifts::Nominal();
576  std::vector<SystShifts> seedShifts = {};
577 
578  //auto bestChi_mtplSeed = fit.Fit(calc, auxShifts, seedPts, seedShifts)->EvalMetricVal();
579  auto bestChi_mtplSeed = fit.Fit(calc, auxShifts, seedPts, seedShifts, IFitter::kQuiet)->EvalMetricVal();
580 
581  if (bestChi_mtplSeed < bestChi) {
582  bestChi = bestChi_mtplSeed;
583  bestDelta13 = kFitDelta13InPiUnitsSterile.GetValue(calc);
584  bestssTh13 = kFitSinSqTheta13Sterile.GetValue(calc);
585  bestssTh23 = kFitSinSqTheta23Sterile.GetValue(calc);
586  bestDmSq = kFitDmSq32Sterile.GetValue(calc);
587  bestDelta24 = kFitDelta24InPiUnitsSterile.GetValue(calc);
588  bestTh24 = kFitTheta24InDegreesSterile.GetValue(calc);
589  bestTh34 = kFitTheta34InDegreesSterile.GetValue(calc);
590  bestDm41 = kFitDmSq41Sterile.GetValue(calc);
591  }
592 
593  SetAngles(calc);
594 
595  }
596  }
597  }
598 
599  else if(dm41vsth24 || dm41vsth34) {
600  for (double th_seed: {1.0e-5, 1.0e-4, 1.0e-3, 1.0e-2, 1.0e-1}) {
601  for (double dm_seed: {5.0e-2, 5.0e-1, 5.0, 5.0e+1}) {
602  std::vector<double> th_seed_vec; th_seed_vec.push_back(th_seed);
603  std::vector<double> dm_seed_vec; dm_seed_vec.push_back(dm_seed);
604  std::map<const IFitVar*, std::vector<double> > seedPts;
605  if(dm41vsth24) {
606  seedPts[&kFitSinSqTheta24Sterile] = th_seed_vec;
607  seedPts[&kFitDmSq41Sterile] = dm_seed_vec;
608  seedPts[&kFitSinSqTheta34Sterile] = {0.5};
609  }
610  else if(dm41vsth34) {
611  seedPts[&kFitSinSqTheta34Sterile] = th_seed_vec;
612  seedPts[&kFitDmSq41Sterile] = dm_seed_vec;
613  seedPts[&kFitSinSqTheta24Sterile] = {0.5};
614  }
615 
616  SystShifts auxShifts = SystShifts::Nominal();
617  std::vector<SystShifts> seedShifts = {};
618 
619  //auto bestChi_mtplSeed = fit.Fit(calc, auxShifts, seedPts, seedShifts)->EvalMetricVal();
620  auto bestChi_mtplSeed = fit.Fit(calc, auxShifts, seedPts, seedShifts, IFitter::kQuiet)->EvalMetricVal();
621 
622  if (bestChi_mtplSeed < bestChi) {
623  bestChi = bestChi_mtplSeed;
624  bestDelta13 = kFitDelta13InPiUnitsSterile.GetValue(calc);
625  bestssTh13 = kFitSinSqTheta13Sterile.GetValue(calc);
626  bestssTh23 = kFitSinSqTheta23Sterile.GetValue(calc);
627  bestDmSq = kFitDmSq32Sterile.GetValue(calc);
628  bestDelta24 = kFitDelta24InPiUnitsSterile.GetValue(calc);
629  bestTh24 = kFitSinSqTheta24Sterile.GetValue(calc);
630  bestTh34 = kFitSinSqTheta34Sterile.GetValue(calc);
631  bestDm41 = kFitDmSq41Sterile.GetValue(calc);
632  }
633 
634  SetAngles(calc);
635 
636  }
637  }
638  }
639 
640  std::clock_t tiend = std::clock();
641  double calctime = ( tiend - tistart ) / (double) CLOCKS_PER_SEC;
642 
643  FCPoint pt(profDelta13_PiUnits, profssTh23, profssTh13, profDmSq,
644  profDelta24_PiUnits, profTh24, profTh34, profDm41,
645  0., 0.,
646  trueDelta13, truessTh23, truessTh13, trueDmSq,
647  trueDelta24, trueTh24, trueTh34, trueDm41,
648  trueChi,
649  bestDelta13, bestssTh23, bestssTh13, bestDmSq,
650  bestDelta24, bestTh24, bestTh34, bestDm41,
651  bestChi,
652  trueChi-bestChi, loadtime, calctime);
653 
654  double bestX = 0;
655  double bestY = 0;
656 
657  if(th24 || th34 || th24vsth34) {
658  bestY = bestTh24;
659  bestX = bestTh34;
660  }
661  else if(dm41vsth24) {
662  bestY = bestDm41;
663  bestX = bestTh24;
664  }
665  else if(dm41vsth34) {
666  bestY = bestDm41;
667  bestX = bestTh34;
668  }
669 
670  std::cout << "Trial: " << iTrial << ", TrueX: " << trueX << ", Y: " << trueY << ", trueChi: " << trueChi << ", bestX: " << bestX << ", Y: " << bestY << ", bestChi: " << bestChi << ", true-best: " << trueChi-bestChi << ", loadTime: " << loadtime << ", calcTime: " << calctime << std::endl;
671 
672  fccol.AddPoint(pt);
673 
674  delete expt;
675  delete mock_spec;
676  }
677  }
678  fccol.SaveToFile(fcout);
679 }
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
void FakeData()
Definition: rootlogon.C:156
const FitSinSqTheta24Sterile kFitSinSqTheta24Sterile
IPrediction * LoadExtrapPrediction(TString opt)
Spectrum LoadCosmicSpectrum(bool rhc=false)
Function to return cosmic spectrum.
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
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
double GetValue(const osc::IOscCalcAdjustable *osc) const override
virtual Spectrum Predict(osc::IOscCalc *calc) const =0
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
void OverridePOT(double newpot)
DO NOT USE UNLESS YOU ARE 110% CERTAIN THERE ISN&#39;T A BETTER WAY!
Definition: Spectrum.h:225
const char * p
Definition: xmltok.h:285
constexpr T pow(T x)
Definition: pow.h:75
std::vector< const ISyst * > LoadISysts()
Function to load ISysts.
void AddPoint(const FCPoint &p)
Definition: FCCollection.h:17
Adapt the PMNS_Sterile calculator to standard interface.
const FitDmSq32Sterile kFitDmSq32Sterile
const double kNus19Th13
Definition: Utilities.h:32
std::vector< double > LogBins(int nbins, double xlo, double xhi)
const FitTheta24InDegreesSterile kFitTheta24InDegreesSterile
const FitSinSqTheta34Sterile kFitSinSqTheta34Sterile
static SystShifts Nominal()
Definition: SystShifts.h:34
osc::OscCalcDumb calc
void SetAngles(osc::OscCalcSterile *calc)
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
void Scale(double c)
Multiply this spectrum by a constant c.
Definition: Spectrum.cxx:260
#define M_PI
Definition: SbMath.h:34
double GetValue(const osc::IOscCalcAdjustable *osc) const override
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const double kAna2018RHCPOT
Definition: Exposures.h:208
const XML_Char const XML_Char * data
Definition: expat.h:268
Represents the results of a single Feldman-Cousins pseudo-experiment.
Definition: FCPoint.h:8
const FitDelta24InPiUnitsSterile kFitDelta24InPiUnitsSterile
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
const int nbins
Definition: cellShifts.C:15
osc::OscCalcSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
Definition: Calcs.cxx:97
expt
Definition: demo5.py:34
const FitSinSqTheta13Sterile kFitSinSqTheta13Sterile
Sum up livetimes from individual cosmic triggers.
std::vector< const IExperiment * > GetConstraints()
Gaussian constrains on atmospheric parameters.
Definition: Utilities.h:416
double GetValue(const osc::IOscCalcAdjustable *osc) const override
virtual std::unique_ptr< IFitSummary > Fit(osc::IOscCalcAdjustable *seed, SystShifts &bestSysts=junkShifts, const SeedList &seedPts=SeedList(), const std::vector< SystShifts > &systSeedPts={}, Verbosity verb=kVerbose) const
Master fitting method. Depends on FitHelper and FitHelperSeeded.
Definition: IFitter.cxx:69
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
std::vector< bool > GetOptionConfig(TString opt, std::vector< std::string > options)
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
Compare a single data spectrum to the MC + cosmics expectation.
Float_t d
Definition: plot.C:236
double GetValue(const osc::IOscCalcAdjustable *osc) const override
void OverrideLivetime(double newlive)
DO NOT USE UNLESS YOU ARE 110% CERTAIN THERE ISN&#39;T A BETTER WAY!
Definition: Spectrum.h:228
const double j
Definition: BetheBloch.cxx:29
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
double GetValue(const osc::IOscCalcAdjustable *osc) const override
const FitDelta13InPiUnitsSterile kFitDelta13InPiUnitsSterile
std::vector< CovarianceMatrix * > GetMatrices(std::vector< covmx::Sample > samples, bool cvmfs=true)
Definition: Utilities.h:350
double GetValue(const osc::IOscCalcAdjustable *osc) const override
std::vector< float > Spectrum
Definition: Constants.h:570
double GetValue(const osc::IOscCalcAdjustable *osc) const override
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
OStream cout
Definition: OStream.cxx:6
std::vector< double > LinBins(int nbins, double xlo, double xhi)
const std::string path
Definition: plot_BEN.C:43
const Binning bins
Definition: NumuCC_CPiBin.h:8
Combine multiple component experiments.
const double kNus19Th23
Definition: Utilities.h:33
::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:323
void make_fc_nus_surfs_nersc_2019(int NPts=1, int startidx=0, int endidx=0, bool nh=true, int N=1, int jid=1, int npts_flush=1, std::string plot="th24vsth34_UO", TString beam="neardetfardetFHCcovmx", bool corrSysts=true, bool readFromFile=true, double max_time_in_sec=3600)
double GetValue(const osc::IOscCalcAdjustable *osc) const override
T sin(T number)
Definition: d0nt_math.hpp:132
IPrediction * LoadPrediction(std::string detector, bool rhc=false, std::string syst_type="all")
Function to load prediction object.
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
const FitSinSqTheta23Sterile kFitSinSqTheta23Sterile
double GetValue(const osc::IOscCalcAdjustable *osc) const override
exit(0)
Base class defining interface for experiments.
Definition: IExperiment.h:14
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
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
double GetValue(const osc::IOscCalcAdjustable *osc) const override
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
const double kNus19Dm21
Definition: Utilities.h:27
const double kAna2018FHCPOT
Definition: Exposures.h:207
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
void SaveToFile(const std::string &fname) const
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
Float_t e
Definition: plot.C:35
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:222
const FitDmSq41Sterile kFitDmSq41Sterile
const double kAna2018FHCLivetime
Definition: Exposures.h:213
const double kNus19SensLivetime
Definition: Utilities.h:39
const double kAna2018RHCLivetime
Definition: Exposures.h:214
std::vector< const ISyst * > GetSystsFromFile(covmx::Sample sample)
Get systematics for a given sample.
Compare a single data spectrum to the MC + cosmics expectation.
void plot(std::string label, std::map< std::string, std::map< std::string, Spectrum * >> &plots, bool log)
const double kNus19Delta13
Definition: Utilities.h:36
Collection of FCPoint. Serializable to/from a file.
Definition: FCCollection.h:12
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17