make_fc_nus_surfs_nersc_2018.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 "CAFAna/Cuts/Cuts.h"
10 #include "CAFAna/FC/FCPoint.h"
11 #include "CAFAna/FC/FCSurface.h"
12 #include "CAFAna/Vars/FitVars.h"
14 #include "CAFAna/Vars/Vars.h"
22 #include "CAFAna/Core/Utilities.h"
23 #include "CAFAna/Core/Progress.h"
24 #include "CAFAna/Core/Loaders.h"
26 #include "CAFAna/Core/IFitVar.h"
27 #include "CAFAna/Analysis/Calcs.h"
29 #include "CAFAna/Fit/Fit.h"
30 #include "CAFAna/Analysis/Style.h"
31 #include "CAFAna/Analysis/Plots.h"
35 
38 
41 
42 #include <sstream>
43 
44 using namespace ana;
45 
46 std::vector<double> LogBins(int nbins,double xlo, double xhi)
47 {
48  double log_spacing = exp( (log(xhi) - log(xlo))/(nbins) );
49  std::vector<double> binning;
50  for (int i = 0; i <= nbins; ++i) {
51  if (i == 0) binning.push_back(xlo);
52  else binning.push_back(log_spacing*binning[i-1]);
53  }
54  return binning;
55 }
56 
57 std::vector<double> LinBins(int nbins,double xlo, double xhi)
58 {
59  double lin_spacing = (xhi - xlo)/(nbins);
60  std::vector<double> binning;
61  for (int i = 0; i <= nbins; ++i) {
62  if (i == 0) binning.push_back(xlo);
63  else binning.push_back(lin_spacing+binning[i-1]);
64  }
65  return binning;
66 }
67 
68 void make_fc_nus_surfs_nersc_2018(int NPts, int startidx, int endidx, bool nh,
69  int N, int jid, int npts_flush = 0,
70  std::string plot = "th24vsth34_UO" ,
71  TString beam = "RHC", bool corrSysts = true,
72  bool readFromFile = false)
73 {
74  //clock starts
75  std::clock_t tstart = std::clock();
76 
77  // Get running options
78  bool fhc = beam.Contains("FHC", TString::kIgnoreCase);
79  bool rhc = beam.Contains("RHC", TString::kIgnoreCase);
80  assert(fhc || rhc);
81 
82  bool covmx = beam.Contains("covmx", TString::kIgnoreCase);
83  bool extrap = beam.Contains("extrap", TString::kIgnoreCase);
84  assert(covmx xor extrap);
85 
86  if (covmx && rhc)
87  throw std::runtime_error("RHC covariance matrix not supported yet.");
88 
89  if (extrap && fhc)
90  throw std::runtime_error("Extrapolated fit not supported for FHC.");
91 
92  std::cout << "beam:: " << beam << std::endl;
93  std::cout << "systs:: " << corrSysts << std::endl;
94  if (covmx && corrSysts) {
95  corrSysts = false;
96  std::cout << "Using covariance matrix, so disabling traditional systematics." << std::endl;
97  }
98 
99  bool upperOctant = false;
100  if(plot == "th24vsth34_UO" ||
101  plot == "th24_UO" ||
102  plot == "th34_UO" ||
103  plot == "dm41vsth24_UO" ||
104  plot == "dm41vsth34_UO")
105  {
106  upperOctant = true;
107  }
108 
109  bool th24 = false;
110  bool th34 = false;
111  bool th24vsth34 = false;
112  bool dm41vsth24 = false;
113  bool dm41vsth34 = false;
114  if(plot == "th24_UO" || plot == "th24_LO"){
115  th24 = true;
116  }
117  else if(plot == "th34_UO" || plot == "th34_LO"){
118  th34 = true;
119  }
120  else if(plot == "th24vsth34_UO" || plot == "th24vsth34_LO") th24vsth34 = true;
121  else if(plot == "dm41vsth24_UO" || plot == "dm41vsth24_LO") dm41vsth24 = true;
122  else if(plot == "dm41vsth34_UO" || plot == "dm41vsth34_LO") dm41vsth34 = true;
123 
124  // bool use_nd;
125  bool use_fd;
126  std::string detectors_in;
127  if(beam.Contains("nd", TString::kIgnoreCase)) {
128  // use_nd = true;
129  use_fd = false;
130  detectors_in = "nd";
131  }
132  else if (beam.Contains("fd", TString::kIgnoreCase) || extrap) {
133  // use_nd = false;
134  use_fd = true;
135  detectors_in = "fd";
136  }
137  else if (beam.Contains("both", TString::kIgnoreCase)) {
138  // use_nd = true;
139  use_fd = true;
140  detectors_in = "both";
141  }
142  else throw std::runtime_error("Detector must be \"nd\", \"fd\" or \"both\".");
143 
144  beam += "extrap";
145  if(corrSysts) beam += "syst";
146 
147  // Get POT & livetimes
149  double fd_pot = rhc ? kAna2018RHCPOT : kAna2018FHCPOT;
150  double fd_livetime = rhc ? kAna2018RHCLivetime : kAna2018FHCLivetime;
151 
152  double pot = use_fd ? fd_pot : nd_pot;
153  double livetime = use_fd ? fd_livetime : 0;
154 
155  IPrediction* pred = nullptr;
156  Spectrum* cosmic = nullptr;
157 
158  if (covmx) {
159 
160  //Default Sim for CovMx Fit
161  pred = GetDefaultSimulation(0, "both", "none");
162  //PredictionConcat* sim = GetDefaultSimulation(rhc, detectors_in); //only fhc "both" available at present
163 
164 
165  //Merged cosmic spectrum
166  std::vector<TH1D*> cosmics = GetCosmics(0, "both");
167  //std::vector<TH1D*> cosmics = GetCosmics(rhc, detectors_in); //only fhc "both" available at present
168  cosmic = new Spectrum(dynamic_cast<PredictionConcat*>(pred)->MergeHists(cosmics), pot, livetime);
169 
170  } // if covmx
171 
172  else {
173 
174  // load prediction for Extrap fit
175  pred = LoadPrediction(beam);
176 
177  // load cosmics for Extrap fit
178  cosmic = new Spectrum(LoadCosmicSpectrum(rhc));
179  cosmic->Scale(livetime/cosmic->Livetime());
180  cosmic->OverridePOT(pot);
181 
182  }
183 
184  // Get the systs
185  std::vector<const ISyst*> systs = {};
186  if (corrSysts)
187  systs = rhc ? GetNus18RHCExtrapSysts() : GetNus18FHCFDSysts("all");
188 
189  // OscCalc for Sterile 3+1 Model
191  calc4f->SetNFlavors(4);
192  SetAngles(calc4f);
193 
194  //Load the covariance matrix generator, make experiment and Gaussian constraints
195  CovarianceMatrix* mx = nullptr;
196  if (covmx) mx = GetFullCovMx(dynamic_cast<PredictionConcat*>(pred), calc4f, pot, "both");
197 
198  // Binnings
199  std::vector<double> th24Bins;
200  std::vector<double> th34Bins;
201  std::vector<double> dm2Bins;
202 
203  if(th24){
204  th24Bins = LinBins(180, 0., 45.);
205  }
206  else if(th34) {
207  th34Bins = LinBins(360, 0., 90.);
208  }
209  else if (th24vsth34) {
210  th24Bins = LinBins(45, 0., 45.);
211  th34Bins = LinBins(90, 0., 90.);
212  }
213  else if (dm41vsth24) {
214  th24Bins = LogBins(50, 1.0e-6, 1);
215  dm2Bins = LogBins(50, 1.0e-2, 100);
216  }
217  else if (dm41vsth34) {
218  th34Bins = LogBins(50, 1.0e-6, 1);
219  dm2Bins = LogBins(50, 1.0e-2, 100);
220  }
221 
222  double binxlowedge = 0.;
223  double binxuppedge = 0.;
224  double binylowedge = 0.;
225  double binyuppedge = 0.;
226  double binxcenter = 0.;
227  double binycenter = 0.;
228 
229  int ibin = startidx;
230  int fbin = endidx;
231 
232  std::string binname = "_bins_"+std::to_string(ibin) + "_to_" + std::to_string(fbin);
233  std::string number = std::to_string(N);
234  std::string hierStr = nh ? "NH" : "IH";
235  std::string fcout = "FCCol_nus18_"+plot+"_"+hierStr+"_";
236  fcout += fhc ? "FHC":"RHC";
237  fcout += corrSysts ? "_systs":"_stats";
238  fcout += binname+"_v"+number+'_'+std::to_string(jid)+".root";
239 
240  // finished loading
241  std::clock_t tload = std::clock();
242  double loadtime = ( tload - tstart ) / (double) CLOCKS_PER_SEC;
243 
244  FCCollection fccol;
245  TRandom3 rnd(0); // zero seeds randomly
246 
247  for (int i = startidx; i <= endidx; i++) {
248 
249  if(th24){
250  binxlowedge = th34Bins.front();
251  binxuppedge = th34Bins.back();
252  binylowedge = th24Bins[i];
253  binyuppedge = th24Bins[i+1];
254  binxcenter = 0.5*(binxlowedge+binxuppedge);
255  binycenter = 0.5*(binylowedge+binyuppedge);
256  }
257  else if(th34){
258  binxlowedge = th34Bins[i];
259  binxuppedge = th34Bins[i+1];
260  binylowedge = th24Bins.front();
261  binyuppedge = th24Bins.back();
262  binxcenter = 0.5*(binxlowedge+binxuppedge);
263  binycenter = 0.5*(binylowedge+binyuppedge);
264  }
265  else if (th24vsth34) {
266  int x_bins = th34Bins.size() - 1;
267  int x_idx = i % x_bins;
268  int y_idx = i / x_bins;
269  binxlowedge = th34Bins[x_idx];
270  binxuppedge = th34Bins[x_idx+1];
271  binylowedge = th24Bins[y_idx];
272  binyuppedge = th24Bins[y_idx+1];
273  binxcenter = 0.5*(binxlowedge+binxuppedge);
274  binycenter = 0.5*(binylowedge+binyuppedge);
275  }
276  else if (dm41vsth24) {
277  int x_bins = th24Bins.size() - 1;
278  int x_idx = i % x_bins;
279  int y_idx = i / x_bins;
280  binxlowedge = th24Bins[x_idx];
281  binxuppedge = th24Bins[x_idx+1];
282  binylowedge = dm2Bins[y_idx];
283  binyuppedge = dm2Bins[y_idx+1];
284  binxcenter = 0.5*(binxlowedge+binxuppedge);
285  binycenter = 0.5*(binylowedge+binyuppedge);
286  }
287  else if (dm41vsth34) {
288  int x_bins = th34Bins.size() - 1;
289  int x_idx = i % x_bins;
290  int y_idx = i / x_bins;
291  binxlowedge = th34Bins[x_idx];
292  binxuppedge = th34Bins[x_idx+1];
293  binylowedge = dm2Bins[y_idx];
294  binyuppedge = dm2Bins[y_idx+1];
295  binxcenter = 0.5*(binxlowedge+binxuppedge);
296  binycenter = 0.5*(binylowedge+binyuppedge);
297  }
298 
299  std::cout << "Bin X (low-center-up): " << binxlowedge << "-" << binxcenter << "-" << binxuppedge << std::endl;
300  std::cout << "Bin Y (low-center-up): " << binylowedge << "-" << binycenter << "-" << binyuppedge << std::endl;
301 
302  double profDelta13_PiUnits = 1.46;
303  double profssTh13 = pow( sin( kNus18Th13 ) , 2);
304  double profssTh23 = 0.470;
305  double profDmSq = kNus18Dm21 + kNus18Dm32;
306  double profDelta24_PiUnits = 1.;
307  if(upperOctant){
308  profDelta13_PiUnits = kNus18Delta13/M_PI;
309  profssTh23 = pow( sin( kNus18Th23 ) , 2);
310  }
311  double profTh24 = 0.;
312  double profTh34 = 0.;
313  double profDm41 = 0.;
314  if(th24 || th34 || th24vsth34) {
315  profTh24 = binxcenter;
316  profTh34 = binycenter;
317  profDm41 = 0.5;
318  }
319  else if(dm41vsth24) {
320  profTh24 = binxcenter;
321  profTh34 = 0.5;
322  profDm41 = binycenter;
323  }
324  else if(dm41vsth34) {
325  profTh24 = 0.5;
326  profTh34 = binxcenter;
327  profDm41 = binycenter;
328  }
329 
330  // Vector of parameters to fit: #theta_{23}, #theta_{34}, #theta_{24}, #deltam_{32}
331  std::vector<const IFitVar*> prof_vars;
332  prof_vars.push_back(&kFitSinSqTheta23Sterile); // to be constrained to PDG
333  prof_vars.push_back(&kFitDmSq32Sterile); // to be constrained to PDG
334  prof_vars.push_back(&kFitDelta24InPiUnitsSterile);
335 
336  if(th24) {
337  prof_vars.push_back(&kFitTheta34InDegreesSterile);
338  }
339  else if(th34) {
340  prof_vars.push_back(&kFitTheta24InDegreesSterile);
341  }
342  else if(dm41vsth24) {
343  prof_vars.push_back(&kFitSinSqTheta34Sterile);
344  }
345  else if(dm41vsth34) {
346  prof_vars.push_back(&kFitSinSqTheta24Sterile);
347  }
348 
349  std::vector<const IFitVar*> fitvars = prof_vars;
350 
351  if(th24) {
352  fitvars.push_back(&kFitTheta24InDegreesSterile);
353  }
354  else if(th34) {
355  fitvars.push_back(&kFitTheta34InDegreesSterile);
356  }
357  else if(th24vsth34) {
358  fitvars.push_back(&kFitTheta24InDegreesSterile);
359  fitvars.push_back(&kFitTheta34InDegreesSterile);
360  }
361  else if(dm41vsth24) {
362  fitvars.push_back(&kFitSinSqTheta24Sterile);
363  fitvars.push_back(&kFitDmSq41Sterile);
364  }
365  else if(dm41vsth34) {
366  fitvars.push_back(&kFitSinSqTheta34Sterile);
367  fitvars.push_back(&kFitDmSq41Sterile);
368  }
369 
370  for(int iTrial = 0; iTrial < NPts; ++iTrial){
371  std::clock_t tistart = std::clock();
372  SetAngles(calc4f);
373 
374  std::cout << "On trial: " << iTrial << std::endl;
375 
376  auto trueX = rnd.Uniform(binxlowedge, binxuppedge);
377  auto trueY = rnd.Uniform(binylowedge, binyuppedge);
378  auto trueDelta24 = rnd.Uniform(0.,2.);
379 
380  double trueTh24 = 0.;
381  double trueTh34 = 0.;
382  double trueDm41 = 0.;
383 
384  int mock_seed = 0;
385  std::cout << "Mock Seed: " << mock_seed << std::endl;
386 
387  if(th24 || th34 || th24vsth34) {
388  trueTh24 = trueY;
389  trueTh34 = trueX;
390  kFitTheta34InDegreesSterile.SetValue(calc4f, trueTh34);
391  kFitTheta24InDegreesSterile.SetValue(calc4f, trueTh24);
392  kFitDelta24InPiUnitsSterile.SetValue(calc4f, trueDelta24);
393  }
394  else if(dm41vsth24) {
395  trueDm41 = trueY;
396  trueTh24 = trueX;
397  trueTh34 = rnd.Uniform(0.,1.);
398  kFitSinSqTheta24Sterile.SetValue(calc4f, trueTh24);
399  kFitDmSq41Sterile.SetValue(calc4f, trueDm41);
400  kFitSinSqTheta34Sterile.SetValue(calc4f, trueTh34);
401  kFitDelta24InPiUnitsSterile.SetValue(calc4f, trueDelta24);
402  }
403  else if(dm41vsth34) {
404  trueDm41 = trueY;
405  trueTh34 = trueX;
406  trueTh24 = rnd.Uniform(0.,1.);
407  kFitSinSqTheta24Sterile.SetValue(calc4f, trueTh24);
408  kFitDmSq41Sterile.SetValue(calc4f, trueDm41);
409  kFitSinSqTheta34Sterile.SetValue(calc4f, trueTh34);
410  kFitDelta24InPiUnitsSterile.SetValue(calc4f, trueDelta24);
411  }
412 
413  // Generate Mock Data
414  Spectrum* mock_spec;
415 
416  // For CovMx Fit
417  if (covmx) {
418 
419  Spectrum s_systs(mx->GetUniverse(pred->Predict(calc4f).ToTH1(pot)), pot, livetime);
420  s_systs += *cosmic;
421  TH1D* hMockSpectrumCovMx = s_systs.MockData(pot,mock_seed).ToTH1(pot);
422  mock_spec = new Spectrum(hMockSpectrumCovMx, pot, livetime);
423  delete hMockSpectrumCovMx;
424 
425  } // if covmx
426 
427  // For Extrap Fit
428  else {
429 
430  Spectrum obs = pred->Predict(calc4f);
431  obs += *cosmic;
432  TH1 *hMockSpectrum = obs.MockData(pot,mock_seed).ToTH1(pot);
433  mock_spec = new Spectrum(hMockSpectrum, pot, livetime);
434  delete hMockSpectrum;
435 
436  } // if extrap
437 
438  // Compares (mock)data sample with MC sample
439  Double_t meanSinSqTheta23Sterile = 0;
440  Double_t sigmaSinSqTheta23Sterile = 0;
441  Double_t meanDmSq32Sterile = 0;
442  Double_t sigmaDmSq32Sterile = 0;
443 
444  meanSinSqTheta23Sterile = profssTh23;
445  sigmaSinSqTheta23Sterile = 0.040;
446 
447  meanDmSq32Sterile = kNus18Dm32;
448  sigmaDmSq32Sterile = kNus18Dm32Sigma;
449 
450  std::vector<const IChiSqExperiment*> expts;
451 
452  GaussianConstraint s2th23Constraint(&kFitSinSqTheta23Sterile, meanSinSqTheta23Sterile, sigmaSinSqTheta23Sterile);
453  GaussianConstraint dmsq32Constraint(&kFitDmSq32Sterile, meanDmSq32Sterile, sigmaDmSq32Sterile);
454 
455  expts.push_back(&s2th23Constraint);
456  expts.push_back(&dmsq32Constraint);
457 
458  const IChiSqExperiment* exp;
459  if (covmx) exp = new OscCovMxExperiment(pred, mx, *mock_spec, *cosmic);
460  else exp = new SingleSampleExperiment(pred, *mock_spec, *cosmic);
461  expts.push_back(exp);
462 
463  MultiExperiment expt(expts);
464 
465  // -------- TRUE FIT -------- //
466  Fitter fittrue(&expt, prof_vars, systs);
467 
468  //double trueChi = fittrue.Fit(calc4f);
469  double trueChi = fittrue.Fit(calc4f, Fitter::kQuiet);
470  double trueDelta13 = kFitDelta13InPiUnitsSterile.GetValue(calc4f);
471  double truessTh13 = kFitSinSqTheta13Sterile.GetValue(calc4f);
472  double truessTh23 = kFitSinSqTheta23Sterile.GetValue(calc4f);
473  double trueDmSq = kFitDmSq32Sterile.GetValue(calc4f);
474  trueDelta24 = kFitDelta24InPiUnitsSterile.GetValue(calc4f);
475 
476  if(th24) {
477  trueTh34 = kFitTheta34InDegreesSterile.GetValue(calc4f);
478  }
479  else if(th34) {
480  trueTh24 = kFitTheta24InDegreesSterile.GetValue(calc4f);
481  }
482  else if(dm41vsth24) {
483  trueTh34 = kFitSinSqTheta34Sterile.GetValue(calc4f);
484  }
485  else if(dm41vsth34) {
486  trueTh24 = kFitSinSqTheta24Sterile.GetValue(calc4f);
487  }
488 
489  // -------- BEST FIT -------- //
490  Fitter fit(&expt, fitvars, systs);
491 
492  double bestChi = 1E20;
493  double bestDelta13 = 0;
494  double bestssTh13 = 0;
495  double bestssTh23 = 0;
496  double bestDmSq = 0;
497  double bestDelta24 = 0;
498  double bestTh24 = 0;
499  double bestTh34 = 0;
500  double bestDm41 = 0;
501 
502  if(th24 || th34 || th24vsth34) {
503  for (double th34_seed: {15, 30, 45, 60, 75}) {
504  for (double th24_seed: {9, 18, 27, 36}) {
505  std::vector<double> th34_seed_vec; th34_seed_vec.push_back(th34_seed);
506  std::vector<double> th24_seed_vec; th24_seed_vec.push_back(th24_seed);
507  std::map<const IFitVar*, std::vector<double> > seedPts;
508  seedPts[&kFitTheta34InDegreesSterile] = th34_seed_vec;
509  seedPts[&kFitTheta24InDegreesSterile] = th24_seed_vec;
510 
511  SystShifts auxShifts = SystShifts::Nominal();
512  std::vector<SystShifts> seedShifts = {};
513 
514  auto bestChi_mtplSeed = fit.Fit(calc4f, auxShifts, seedPts, seedShifts);
515  // auto bestChi_mtplSeed = fit.Fit(calc4f, auxShifts, seedPts, seedShifts, Fitter::kQuiet);
516 
517  if (bestChi_mtplSeed < bestChi) {
518  bestChi = bestChi_mtplSeed;
519  bestDelta13 = kFitDelta13InPiUnitsSterile.GetValue(calc4f);
520  bestssTh13 = kFitSinSqTheta13Sterile.GetValue(calc4f);
521  bestssTh23 = kFitSinSqTheta23Sterile.GetValue(calc4f);
522  bestDmSq = kFitDmSq32Sterile.GetValue(calc4f);
523  bestDelta24 = kFitDelta24InPiUnitsSterile.GetValue(calc4f);
524  bestTh24 = kFitTheta24InDegreesSterile.GetValue(calc4f);
525  bestTh34 = kFitTheta34InDegreesSterile.GetValue(calc4f);
526  bestDm41 = 0.5;
527  }
528 
529  SetAngles(calc4f);
530 
531  }
532  }
533  }
534 
535  else if(dm41vsth24 || dm41vsth34) {
536  for (double th_seed: {1.0e-5, 1.0e-4, 1.0e-3, 1.0e-2, 1.0e-1}) {
537  for (double dm_seed: {5.0e-2, 5.0e-1, 5.0, 5.0e+1}) {
538  std::vector<double> th_seed_vec; th_seed_vec.push_back(th_seed);
539  std::vector<double> dm_seed_vec; dm_seed_vec.push_back(dm_seed);
540  std::map<const IFitVar*, std::vector<double> > seedPts;
541  if(dm41vsth24) {
542  seedPts[&kFitSinSqTheta24Sterile] = th_seed_vec;
543  seedPts[&kFitDmSq41Sterile] = dm_seed_vec;
544  seedPts[&kFitSinSqTheta34Sterile] = {0.5};
545  }
546  else if(dm41vsth34) {
547  seedPts[&kFitSinSqTheta34Sterile] = th_seed_vec;
548  seedPts[&kFitDmSq41Sterile] = dm_seed_vec;
549  seedPts[&kFitSinSqTheta24Sterile] = {0.5};
550  }
551 
552  SystShifts auxShifts = SystShifts::Nominal();
553  std::vector<SystShifts> seedShifts = {};
554 
555  auto bestChi_mtplSeed = fit.Fit(calc4f, auxShifts, seedPts, seedShifts);
556  // auto bestChi_mtplSeed = fit.Fit(calc4f, auxShifts, seedPts, seedShifts, Fitter::kQuiet);
557 
558  if (bestChi_mtplSeed < bestChi) {
559  bestChi = bestChi_mtplSeed;
560  bestDelta13 = kFitDelta13InPiUnitsSterile.GetValue(calc4f);
561  bestssTh13 = kFitSinSqTheta13Sterile.GetValue(calc4f);
562  bestssTh23 = kFitSinSqTheta23Sterile.GetValue(calc4f);
563  bestDmSq = kFitDmSq32Sterile.GetValue(calc4f);
564  bestDelta24 = kFitDelta24InPiUnitsSterile.GetValue(calc4f);
565  bestTh24 = kFitSinSqTheta24Sterile.GetValue(calc4f);
566  bestTh34 = kFitSinSqTheta34Sterile.GetValue(calc4f);
567  bestDm41 = kFitDmSq41Sterile.GetValue(calc4f);
568  }
569 
570  SetAngles(calc4f);
571 
572  }
573  }
574  }
575 
576  std::clock_t tiend = std::clock();
577  double calctime = ( tiend - tistart ) / (double) CLOCKS_PER_SEC;
578 
579  FCPoint pt(profDelta13_PiUnits, profssTh23, profssTh13, profDmSq,
580  profDelta24_PiUnits, profTh34, profTh24, profDm41,
581  0., 0.,
582  trueDelta13, truessTh23, truessTh13, trueDmSq,
583  trueDelta24, trueTh34, trueTh24, trueDm41,
584  trueChi,
585  bestDelta13, bestssTh23, bestssTh13, bestDmSq,
586  bestDelta24, bestTh34, bestTh24, 0.5,
587  bestChi,
588  trueChi-bestChi, loadtime, calctime);
589 
590  double bestX = 0;
591  double bestY = 0;
592 
593  if(th24 || th34 || th24vsth34) {
594  bestY = bestTh24;
595  bestX = bestTh34;
596  }
597  else if(dm41vsth24) {
598  bestY = bestDm41;
599  bestX = bestTh24;
600  }
601  else if(dm41vsth34) {
602  bestY = bestDm41;
603  bestX = bestTh34;
604  }
605 
606  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;
607 
608  fccol.AddPoint(pt);
609  unsigned int nfc_pts = fccol.NPoints();
610 
611  if(npts_flush != 0) {
612 
613  if(nfc_pts % npts_flush == 0 ) {
614  std::cout << "Saving " + std::to_string(nfc_pts) + "th point to " << fcout << std::endl;
615  fccol.SaveToFile(fcout);
616  }
617  }
618 
619  delete exp;
620  }
621  }
622 
623  fccol.SaveToFile(fcout);
624 
625 }
std::vector< const ISyst * > GetNus18FHCFDSysts(std::string syst_type)
Definition: Nus18Systs.cxx:630
const FitSinSqTheta24Sterile kFitSinSqTheta24Sterile
TH1D * GetUniverse(const TH1D *nominal)
Spectrum LoadCosmicSpectrum(bool rhc=false)
Function to return cosmic spectrum.
Oscillation analysis framework, runs over CAF files outside of ART.
double GetValue(const osc::IOscCalculatorAdjustable *osc) const override
std::vector< SystGroupDef > systs
Definition: syst_header.h:384
Adapt the PMNS_Sterile calculator to standard interface.
void SetValue(osc::IOscCalculatorAdjustable *osc, double val) const override
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:553
double GetValue(const osc::IOscCalculatorAdjustable *osc) const override
A simple Gaussian constraint on an arbitrary IFitVar.
const double kNus18Delta13
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:269
constexpr T pow(T x)
Definition: pow.h:75
const double kNus18Dm21
std::vector< double > LinBins(int nbins, double xlo, double xhi)
void AddPoint(const FCPoint &p)
Definition: FCCollection.h:17
const FitDmSq32Sterile kFitDmSq32Sterile
void plot()
double GetValue(const osc::IOscCalculatorAdjustable *osc) const override
const double kAna2018SensitivityFHCNDPOT
Definition: Exposures.h:210
const FitTheta24InDegreesSterile kFitTheta24InDegreesSterile
const FitSinSqTheta34Sterile kFitSinSqTheta34Sterile
static SystShifts Nominal()
Definition: SystShifts.h:34
const double kNus18Th13
double GetValue(const osc::IOscCalculatorAdjustable *osc) const override
void Scale(double c)
Multiply this spectrum by a constant c.
Definition: Spectrum.cxx:734
#define M_PI
Definition: SbMath.h:34
double GetValue(const osc::IOscCalculatorAdjustable *osc) const override
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
unsigned int NPoints() const
Definition: FCCollection.h:26
double GetValue(const osc::IOscCalculatorAdjustable *osc) const override
double GetValue(const osc::IOscCalculatorAdjustable *osc) const override
Represents the results of a single Feldman-Cousins pseudo-experiment.
Definition: FCPoint.h:8
const FitDelta24InPiUnitsSterile kFitDelta24InPiUnitsSterile
const int nbins
Definition: cellShifts.C:15
expt
Definition: demo5.py:34
const FitSinSqTheta13Sterile kFitSinSqTheta13Sterile
Sum up livetimes from individual cosmic triggers.
void SetValue(osc::IOscCalculatorAdjustable *osc, double val) const override
void SetValue(osc::IOscCalculatorAdjustable *osc, double val) const override
void SetAngles(osc::OscCalculatorSterile *calc)
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
Base class defining interface for experiments.
void SetValue(osc::IOscCalculatorAdjustable *osc, double val) const override
#define pot
Compare a single data spectrum to the MC + cosmics expectation.
const double kNus18Dm32
const double kNus18Th23
const FitDelta13InPiUnitsSterile kFitDelta13InPiUnitsSterile
double GetValue(const osc::IOscCalculatorAdjustable *osc) const override
OStream cout
Definition: OStream.cxx:6
A class for generating a covariance matrices as a function of oscillation parameters and systematics ...
Combine multiple component experiments.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void SetValue(osc::IOscCalculatorAdjustable *osc, double val) const override
Spectrum MockData(double pot, int idx=0) const
Mock data is FakeData with Poisson fluctuations applied.
Definition: Spectrum.cxx:787
void SetValue(osc::IOscCalculatorAdjustable *osc, double val) const override
std::vector< const ISyst * > GetNus18RHCExtrapSysts()
Definition: Nus18Systs.cxx:727
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.
const FitSinSqTheta23Sterile kFitSinSqTheta23Sterile
const double kAna2018SensitivityRHCNDPOT
Definition: Exposures.h:211
std::vector< double > LogBins(int nbins, double xlo, double xhi)
std::vector< TH1D * > GetCosmics(bool rhc, std::string detector)
Function to get default cosmics.
const FitTheta34InDegreesSterile kFitTheta34InDegreesSterile
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
const double kAna2018FHCPOT
Definition: Exposures.h:207
osc::OscCalculatorSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
Definition: Calcs.cxx:97
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
void SaveToFile(const std::string &fname) const
Float_t e
Definition: plot.C:35
double GetValue(const osc::IOscCalculatorAdjustable *osc) const override
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:266
const double kNus18Dm32Sigma
const FitDmSq41Sterile kFitDmSq41Sterile
const double kAna2018FHCLivetime
Definition: Exposures.h:213
void make_fc_nus_surfs_nersc_2018(int NPts, int startidx, int endidx, bool nh, int N, int jid, int npts_flush=0, std::string plot="th24vsth34_UO", TString beam="RHC", bool corrSysts=true, bool readFromFile=false)
const double kAna2018RHCLivetime
Definition: Exposures.h:214
double GetValue(const osc::IOscCalculatorAdjustable *osc) const override
Compare a single data spectrum to the MC + cosmics expectation.
Collection of FCPoint. Serializable to/from a file.
Definition: FCCollection.h:12