make_fc_slices_nersc_2018.C
Go to the documentation of this file.
1 #include <fstream>
2 
4 #include "CAFAna/Core/Spectrum.h"
6 
9 
13 #include "CAFAna/Vars/FitVars.h"
14 
15 #include "CAFAna/FC/FCPoint.h"
16 #include "CAFAna/FC/FCCollection.h"
17 #include "CAFAna/FC/FCSurface.h"
18 
22 
24 
25 #include "CAFAna/Analysis/Calcs.h"
26 #include "OscLib/IOscCalc.h"
27 
28 #include "TRandom3.h"
29 #include "TH1.h"
30 #include "TGraph.h"
31 
32 using namespace ana;
33 
34 void SetSeeds(osc::IOscCalcAdjustable* calc, double delta, double th23, double dmsq)
35 {
37  kFitDeltaInPiUnits.SetValue(calc, delta);
38  kFitSinSqTheta23.SetValue(calc, th23);
39  kFitDmSq32.SetValue(calc, dmsq);
40 }
41 
42 // Conventionally, we run with NPts = 100 in each bin for std production of FC points
43 // This script is intended for use on NERSC where we run over a range of bins with a few points each,
44 // sending out multiple instances to get the total 100 points in each bin.
45 // We've usually set 80 bins for slice plots, so bin runs from 0-79.
46 // N is an identifier marker to add to the job output, increment if you ever
47 // need to top up stats anywhere so you don't overwrite same-named files.
48 // plot is typically either delta, dmsq32, or ssth23. Per special request, it is also
49 // able to handle th13 slice but we never show this to the public. It was mainly an
50 // exercise to test NERSC implementation.
51 void make_fc_slices_nersc_2018(int NPts, int startidx, int endidx, bool nh,
52  int N, int jid, int npts_flush = 0,
53  TString plot="delta_LO", TString beam = "both",
54  bool corrSysts = false,
55  bool test = false)
56 {
57  // We want our slices this year to include octant slices
58  // UO: Upper Octant
59  // LO: Lower Octant
60  assert(plot == "delta_UO" || plot == "delta_LO"
61  || plot == "ssth23"
62  /*|| plot == "th13_UO" || plot == "th13_LO" */ // th13 not ready yet
63  || plot == "dmsq_UO" || plot == "dmsq_LO"
64  || plot == "th13_UO" || plot == "th13_LO");
65 
66  assert(beam == "both" || beam == "RHCOnly" || beam == "FHCOnly");
67 
68  bool RHCOnly = beam.Contains("RHCOnly");
69  bool FHCOnly = beam.Contains("FHCOnly");
70  bool both = beam.Contains("both");
71  bool upperOctant = plot.Contains("UO");
72  bool deltaSlice = plot.Contains("delta");
73  bool ssth23Slice = plot.Contains("ssth23");
74  bool th13Slice = plot.Contains("th13");
75  bool dmsqSlice = plot.Contains("dmsq");
76  // Filling the array of bins to FC from previously generated text files.
77  ifstream binlist;
78  std::string hier = nh ? "_NH_": "_IH_";
79 
80  int NBins = 60;
81  int bins[NBins];
82  std::memset(bins,0,NBins);
83 
84  // sequential array of bins to FC Correct
85  for(int i = 0; i < NBins; i++) bins[i] = i;
86 
87  // What bins are we covering in this job?
88  int ibin = bins[startidx];
89  int fbin = bins[endidx];
90 
91  std::string number = std::to_string(N);
92  std::string plot_str;
93  if(ssth23Slice) plot_str = "ssth23";
94  if(deltaSlice) plot_str = "delta";
95  if(dmsqSlice) plot_str = "dmsq";
96  if(th13Slice) plot_str = "th13";
97  if(!ssth23Slice) plot_str += upperOctant? "UO" : "LO";
98  std::string fcout = "FCCol_"+plot_str+hier;
99  fcout += "_"+beam;
100  fcout += corrSysts? "_systs_":"_stats_";
101  fcout +="bins_"+std::to_string(ibin)+"_to_"+std::to_string(fbin)+
102  "_v"+number+'_'+std::to_string(jid)+".root";
103 
104  FCCollection fccol;
106 
107 
108  // We're setting up the requisite inputs to the nue experiment
109 
110  ///////////////////////////////////////////////////////////////////
111  // Nue Prediction
112  ///////////////////////////////////////////////////////////////////
113 
114  double thisPOT;
115  double thisLivetime;
116 
117  // Grab the systs
118  // bools to make the function calls a little more readable
119  std::vector<const ISyst*> systs = {};
120 
121  if(corrSysts) systs = getAllAna2018Systs(kJointAna2018, true, kBoth, true);
122 
123  // Beam mode Norm systs are annoying.... Have to make two systs lists so that
124  // FHC doesn't pick up RHC Normalization syst. Hopefully this will be changed soon
125  // r1548
126  std::vector<const ISyst*> snue_fhc = {};
127  if(corrSysts) snue_fhc = getAllAna2018Systs(kNueAna2018, true, kFHC);
128  std::vector<const ISyst*> snue_rhc = {};
129  if(corrSysts) snue_rhc = getAllAna2018Systs(kNueAna2018, true, kRHC);
130  std::vector<const ISyst*> snumu = {};
131  if(corrSysts) snumu = getAllAna2018Systs(kNumuAna2018);
132  // We're setting up the requisite inputs to the nue experiment
133  // Nue Prediction
134 
135  std::vector <const IPrediction *> nue_preds;
136  std::vector <std::pair <TH1D*, double > > nue_cosmics;
137  std::vector <Spectrum> nue_cosmics_spectra;
138  bool minimizeMemory = true;
139  // Grab the nue predictions for fhc and rhc
140 
141  if(!RHCOnly)
142  // nue_preds.push_back(GetNuePrediction2018("combo", calc, corrSysts, "fhc", false, false, minimizeMemory, true));
143  nue_preds.push_back(GetNuePrediction2018("combo", calc, corrSysts, "fhc", false, false, minimizeMemory, true));
144 
145  if(!FHCOnly)
146  // nue_preds.push_back(GetNuePrediction2018("prop", calc, corrSysts, "rhc", false, false, minimizeMemory, true));
147  nue_preds.push_back(GetNuePrediction2018("prop", calc, corrSysts, "rhc", false, false, minimizeMemory, true));
148 
149  std::cout << "\n------------Done getting Nue Predictions--------------\n" ;
150 
151  if(!RHCOnly) {
152  nue_cosmics.push_back(GetNueCosmics2018("fhc", false, true));
153  nue_cosmics_spectra.push_back(Spectrum(nue_cosmics.back().first,
155  }
156  if(!FHCOnly){
157  nue_cosmics.push_back(GetNueCosmics2018("rhc", false, true));
158  nue_cosmics_spectra.push_back(Spectrum(nue_cosmics.back().first,
160  }
161 
162 
163  //////////////////////////////////////////////////////////////////////
164  // Numu Prediction
165  //////////////////////////////////////////////////////////////////////
166  std::vector<const IPrediction*> numu_preds;
167  std::vector< std::pair<TH1D*, double> > numu_cosmics;
168  std::vector<Spectrum> numu_cosmics_spectra;
169 
170  std::vector<const IPrediction*> this_numu_preds;
171  std::vector< std::pair<TH1D*, double> > this_numu_cosmics;
172 
173  std::vector<std::string> beams;
174  if(!RHCOnly) beams.push_back("fhc");
175  if(!FHCOnly) beams.push_back("rhc");
176  // Construct numu inputs
177 
178  std::cout << "Loading Numu predictions for " << beam << endl;
179  for(std::string beam : beams){
180  this_numu_preds = GetNumuPredictions2018(4,corrSysts, beam, false, kNuMu, minimizeMemory, true);
181  this_numu_cosmics = GetNumuCosmics2018(4, beam, false, true);
182  numu_preds.insert(numu_preds.end(),
183  this_numu_preds.begin(),
184  this_numu_preds.end());
185  numu_cosmics.insert(numu_cosmics.end(),
186  this_numu_cosmics.begin(),
187  this_numu_cosmics.end());
188  for(std::pair<TH1D*, double> h : this_numu_cosmics){
189  thisPOT = kAna2018FHCPOT;
190  thisLivetime = kAna2018FHCLivetime;
191  if(beam == "rhc") {
192  thisPOT = kAna2018RHCPOT;
193  thisLivetime = kAna2018RHCLivetime;
194  }
195  numu_cosmics_spectra.push_back(Spectrum(h.first, thisPOT, thisLivetime));
196  }
197 
198  }
199 
200  std::cout << "----------------Done getting numu predictions--------------\n";
201  // Grab analysis slices to find bin center
202  // TODO clean this up
203  TString s = corrSysts? "systs" : "stats";
204  TString helperupsname = std::string(getenv("FCHELPERANA2018_LIB_PATH"));
205  TString anaslicename = helperupsname
206  + "/slices/hist_slices_2017_joint_realData_"+beam+"combo_"+s;
207 
208  if(th13Slice) anaslicename += "_th13";
209  if(ssth23Slice) anaslicename += "_th23_noOct";
210  if(dmsqSlice) anaslicename += "_dmsq";
211 
212  anaslicename += test? "_test" : "";
213  anaslicename += ".root";
214 
215  // slice names changed in FCHELPERANA2018_LIB_PATH_v1.0
216  /* anaslicename += "_joint_both_fake2017";
217  anaslicename += test? "_noOct_systs" : "_systs_noOct"; // HACK
218  anaslicename += ".root"; */
219 
220  TFile *anaslicefile = TFile::Open(anaslicename);
221  std::cout << "Opened anaslice file " << endl;
222  anaslicefile->cd();
223  std::string plotShortName="delta";
224  if(th13Slice) plotShortName = "th13";
225  if(ssth23Slice) plotShortName = "th23";
226  if(dmsqSlice) plotShortName = "dmsq";
227  plotShortName += nh? "_NH" : "_IH";
228  if(!ssth23Slice)
229  plotShortName += upperOctant? "UO" : "LO";
230 
231  TH1* anaslice = (TH1F*)anaslicefile->Get(("slice_"+ plotShortName).c_str());
232  int NX = anaslice->GetXaxis()->GetNbins();
233  double minX = anaslice->GetXaxis()->GetBinLowEdge(1); // root bin starts at 1
234  double maxX = anaslice->GetXaxis()->GetBinUpEdge(NX);
235  double widthX = anaslice->GetXaxis()->GetBinWidth(1);
236  std::cout << "Got ana histogram" << endl;
237  assert( NBins == NX && "Number of bins in bin array does not match slice");
238 
239  // Grab the FCHelper2017 file
240  // This file saved the best fit parameters at every value ofthe variable
241  // we're plotting
242  TString helpername =
243  helperupsname
244  + "/FCInput/"
245  + "/slices_FCInput_2018_";
246  if(deltaSlice) {
247  if(corrSysts) helpername += "combojoint_realData_"+beam+"combo_"+s+".root";
248  else helpername += "joint_realData_"+beam+"combo_stats.root";
249  }
250 
251  if(ssth23Slice) helpername += "joint_realData_"+beam+"combo_"+s+"_th23_noOct.root";
252  if(dmsqSlice) helpername += "joint_realData_"+beam+"combo_"+s+"_dmsq.root";
253  if(th13Slice) helpername += "joint_realData_"+beam+"combo_"+s+"_th13.root";
254 
255  TFile *fchelp = TFile::Open(helpername);
256  std::cout << "Opened FC helper " << endl;
257  // Let's start throwing
258  int bin;
259  for(int i = startidx; i <= endidx; i++){
260  bin = bins[i];
261  std::cout << bin << " ";
262  std::cout<<"- "<<i<<std::endl;
263 
264 
265 
266  // bin array is initialized with zeros before read from binlist
267  // so there will be trailing zeros in the array after the last bin
268  // from the binlist has filled.
269  // Once we run out of important bins, get out.
270  if(bin == 0 && i != 0) {std::cout << "Ran out of important bins. Continuing.\n"; continue;}
271 
272  anaslicefile->cd();
273  double centerX = anaslice->GetXaxis()->GetBinCenter(bin+1);
274  fchelp->cd();
275 
276  // A note on throwing experiments - we've always set the nuisance parameters
277  // for plots, the parameters in the fit that don't show up in the plot, to
278  // the best fit at each given value of the intersting physics parameter
279  // you're plotting (the frequentist way of doing this). To do that, we
280  // need to store the best fit values of all nuisance parameters as a function
281  // of the physics parameter - that's done with TGraph's are a function of the
282  // physics parameter and there's one stored for each hierarchy / nuisance
283  // parameter combo
284 
285 
286  // Need to look up the right hierarchy when searching for best fits
287  std::string hiername = nh ? "NH" : "IH";
288  std::string oct = upperOctant? "UO": "LO";
289  if(ssth23Slice) oct = "";
290  // Pull out hists and get values at current bin
291  // Do this for the oscillation parameters first
292  // Be sure to not find a seed for the variable you're FC'ing!!!
293  double seeddelta = 0;
294  if (!deltaSlice){
295  // HACK HACK HACK: someone changed the name of things on me..
296  std::string tempDeltaInPi = "_DeltaCPpi";
297  if (th13Slice) tempDeltaInPi = "_DeltaCPpi";
298  TGraph *hdelta=(TGraph*)fchelp->Get((hiername+oct+tempDeltaInPi).c_str());
299  seeddelta = hdelta->Eval(centerX);
300  delete hdelta;
301  }
302  double seedss2th13 = 0.085;
303  if (!th13Slice){
304  TGraph *hss2th13=(TGraph*)fchelp->Get((hiername+oct+"_SinSq2Theta13").c_str());
305  seedss2th13 = hss2th13->Eval(centerX);
306  delete hss2th13;
307  }
308  double seedssth23 = 0.6;
309  if (!ssth23Slice){
310  TGraph *hssth23=(TGraph*)fchelp->Get((hiername+oct+"_SinSqTheta23").c_str());
311  seedssth23 = hssth23->Eval(centerX);
312  delete hssth23;
313  }
314  double seeddmsq32 = hiername == "NH" ? 0.00244 : -0.00244;;
315  if (!dmsqSlice){
316  TGraph *hdmsq32=(TGraph*)fchelp->Get((hiername+oct+"_DmSq32").c_str());
317  seeddmsq32 = hdmsq32->Eval(centerX);
318  delete hdmsq32;
319  }
320 
321  // Now, look for all the systematics
322  std::map<std::string,double> seedsyst = {};
323  for (const ISyst* syst : systs){
324  TGraph* h = (TGraph*)fchelp->Get((hiername+oct+"_"+syst->ShortName()).c_str());
325  if (!h){
326  std::cout << "Don't see a prof hist for " << syst->ShortName() << ". Continuing, but check your ups version." << std::endl;
327  continue;
328  }
329  double seedval = h->Eval(centerX);
330  // std::cout << "seeding " << syst->ShortName() << " at " << seedval << std::endl;
331  seedsyst.insert(std::pair<std::string,double>(syst->ShortName(),seedval));
332  }
333 
334 
335  // Need random numbers to sample the bin
336  TRandom3 rnd(0);
337 
338 
339  // Now everything is set up, we can start throwing FC mock experiments
340  for (int i = 0; i < NPts; i++){
341  ResetOscCalcToDefault(calc);
342  std::cout << "Throw " << std::to_string(i) << endl;
343 
344  double trueX = centerX;
345  std::cout << "truex " << trueX << std::endl;
346  // since we're doing octant and hierarchy slices,
347  // we'll need fit vars that are constrained to the octant/hierarchy
348  // that we're correcting
349  const IFitVar* th23_var;
350  if(ssth23Slice) th23_var = &kFitSinSqTheta23;
351  else{
352  if(upperOctant) th23_var = kFitSinSqTheta23UpperOctant;
353  else th23_var = kFitSinSqTheta23LowerOctant;
354  }
355  const IFitVar* dmsq_var;
356  if(nh) dmsq_var = &kFitDmSq32NormalHierarchy;
357  else dmsq_var = &kFitDmSq32InvertedHierarchy;
358 
359  // The variable you're plotting isn't a nuisance parameter, should
360  // be set to the trueX that you threw for this exeriment.
361  // Other oscillation parameters are nuisance parameters and should be
362  // seeded at their best fit
363  kFitDeltaInPiUnits.SetValue(calc, deltaSlice ?trueX:seeddelta);
364  th23_var->SetValue(calc, ssth23Slice ?trueX:seedssth23);
365  kFitSinSq2Theta13. SetValue(calc, th13Slice ?trueX:seedss2th13);
366  dmsq_var->SetValue(calc, dmsqSlice ?trueX:seeddmsq32);
367  // Set syst shift - they're nuisance parameters, so set to the best fit
368  // Nue and numu predictions are taught about slightly different sets of
369  // systematics, so we need to throw the predictions with different seeds
370  SystShifts seednue_fhc;
371  for (const ISyst *syst : snue_fhc){
372  if (seedsyst.find(syst->ShortName()) == seedsyst.end()) continue;
373  seednue_fhc.SetShift(syst,seedsyst[syst->ShortName()]);
374  }
375  SystShifts seednue_rhc;
376  for (const ISyst *syst : snue_rhc){
377  if (seedsyst.find(syst->ShortName()) == seedsyst.end()) continue;
378  seednue_rhc.SetShift(syst,seedsyst[syst->ShortName()]);
379  }
380 
381  SystShifts seednumu;
382  for (const ISyst *syst : snumu){
383  if (seedsyst.find(syst->ShortName()) == seedsyst.end()) continue;
384  seednumu.SetShift(syst,seedsyst[syst->ShortName()]);
385  }
386 
387 
388 
389  // Make a mock data spectrum for nue
390  // Use mock data + cosmic as predictions
391  // assert(nue_preds.size() == 2 && "Something is wrong with nue preds. Size != 2");
392  std::vector<Spectrum> mocknue;
393  SystShifts this_seednue;
394  int mock_integral = 0;
395  double thisPOT;
396  for (int j = 0; j < (int) nue_preds.size(); j++){
397 
398  if(j==0 && !RHCOnly) {
399  this_seednue = seednue_fhc;
400  thisPOT = kAna2018FHCPOT;
401  }
402  else {
403  this_seednue = seednue_rhc;
404  thisPOT = kAna2018RHCPOT;
405  }
406 
407  Spectrum fakenue = nue_preds[j]->Predict(calc);
408  if(corrSysts) fakenue = nue_preds[j]->PredictSyst(calc,this_seednue);
409  // Add in cosmic prediction here, make sure it gets Poisson fluctuated
410  fakenue += nue_cosmics_spectra[j];
411  mocknue.push_back(fakenue.MockData(thisPOT));
412  // Save this, for posterity's sake
413  mock_integral += mocknue[j].Integral(thisPOT);
414  }
415  std::cout << "---------------Made "
416  << nue_preds.size()
417  << " mock nue data spectra-------------\n";
418 
419  std::vector<Spectrum> mocknumu;
420  // Make a mock data spectrum for numu
421  for (int j = 0; j < (int)numu_preds.size(); j++){
422  double thisPOT = 0;
423  if(j < 3 && !RHCOnly){
424  thisPOT = kAna2018FHCPOT;
425  }
426  else {
427  thisPOT = kAna2018RHCPOT;
428  }
429  Spectrum fakenm = numu_preds[j]->Predict(calc);
430  if(corrSysts) fakenm = numu_preds[j]->PredictSyst(calc,seednumu);
431 
432  fakenm += numu_cosmics_spectra[j]; // Again, add cosmic spectra pre-fluctuation
433  mocknumu.push_back(fakenm.MockData(thisPOT)); // Fluctuate
434  }
435  std::cout << "---------------Made "
436  << numu_preds.size()
437  << " mock numu data spectra-------------\n";
438 
439  // Set up experiment, nova nue + nova numu (4 quantiles) + PDG reactor
440  std::vector<const IExperiment*> expts;
441  std::vector<const IExperiment*> numu_expts;
442  // Nue experiments
443  for(int j = 0; j < (int)nue_preds.size(); j++) {
444  expts.push_back(new SingleSampleExperiment(nue_preds[j], mocknue[j],
445  nue_cosmics[j].first,
446  nue_cosmics[j].second)); // Nue
447  }
448  for (int j = 0; j < (int)numu_preds.size(); j++) {
449  expts.push_back(new SingleSampleExperiment(numu_preds[j], mocknumu[j],
450  numu_cosmics[j].first,
451  numu_cosmics[j].second));
452  numu_expts.push_back(new SingleSampleExperiment(numu_preds[j], mocknumu[j],
453  numu_cosmics[j].first,
454  numu_cosmics[j].second));
455  }
456 
457  if (plot != "th13") {
458  expts.push_back( WorldReactorConstraint2017() );
459  numu_expts.push_back( WorldReactorConstraint2017() );
460  }
461  std::cout << "Creating multiexperiment from " << expts.size() << " experiments" << endl;
462  MultiExperiment expt(expts);
463  MultiExperiment numu_expt(numu_expts);
464 
465  // Correlate your systs
466  // A bit tricky since nue and numu predictions know about a different list
467  // of systs, but there are functions to do this in joint_fit_2017_tools.h
468  for(int j = 0; j < (int)nue_preds.size(); j++) {
469  if(j==0 && !RHCOnly)
470  expt.SetSystCorrelations(j, GetCorrelations(true, true)); // nue FHC
471  else
472  expt.SetSystCorrelations(j, GetCorrelations(true, false)); // nue RHC
473  }
474  auto notfornumu = GetCorrelations(false, false); // Numu not dependent
475  // on beam mode
476  for(int j = 0; j < (int)numu_preds.size(); ++j) {
477  expt.SetSystCorrelations(j+nue_preds.size(), notfornumu);
478  }
479 
480  // Finally, we're get to do some FC
481  // Set up fitter for #chi^{2}_{best}
482 
483  // make a vector of oscillation parameters that are our nuisance parameters
484  std::vector<const IFitVar*> trueVars;
485  if (deltaSlice){
486  trueVars.push_back(dmsq_var);
487  trueVars.push_back(&kFitSinSq2Theta13); // no delta included in this case
488  trueVars.push_back(th23_var);
489  }
490  if (ssth23Slice){
491  trueVars.push_back(dmsq_var);
492  trueVars.push_back(&kFitDeltaInPiUnits);
493  trueVars.push_back(&kFitSinSq2Theta13);
494  }
495  if (dmsqSlice){
496  trueVars.push_back(&kFitDeltaInPiUnits);
497  trueVars.push_back(&kFitSinSq2Theta13);
498  trueVars.push_back(th23_var);
499  }
500  if (th13Slice){
501  trueVars.push_back(&kFitDeltaInPiUnits);
502  trueVars.push_back(dmsq_var);
503  trueVars.push_back(th23_var);
504 
505  }
506  ///////////////////////////////////////////////////////////////////
507  // Numu only fit to find the ssth23 and dmsq seeds
508  //////////////////////////////////////////////////////////////////
509  // in case some systs need different seeds
510 
511  std::vector<SystShifts> junkSystShifts = {};
512  std::vector<SystShifts> thisSystShift= {};
513  if(corrSysts && deltaSlice){
514  if(!RHCOnly){ // if fhc mode, add cherenkov
515  for(double systshift : {-.5, .5}){
516  SystShifts tempShifts (&kAnaCherenkovSyst, systshift);
517  thisSystShift.push_back(tempShifts);
518  }
519  }
520  if(both){ // add neutron only if using both beam modes
521  for(double systshift : {-.5, .5}){
522  SystShifts tempShifts (&kNeutronVisEScalePrimariesSyst2018, systshift);
523  thisSystShift.push_back(tempShifts);
524  }
525  }
526 
527  }
528 
529  std::cout << "Fitting numu to find pre seeds" << endl;
530 
532  SystShifts auxShifts = SystShifts::Nominal();
533  double maxmixing = 0.514;
534 
535  MinuitFitter fitnumu_only(&numu_expt, {&kFitSinSqTheta23, &kFitDmSq32});
536 
537  double minchi_numu = fitnumu_only.Fit(numu_calc, auxShifts,
538  {{&kFitDmSq32, {nh? 2.35e-3 : -2.35e-3}},
539  {&kFitSinSqTheta23, {0.4}}},
540  junkSystShifts, IFitter::kQuiet)->EvalMetricVal();
541 
542  double pre_seed_th23 = kFitSinSqTheta23.GetValue(numu_calc);
543  double pre_seed_dmsq = kFitDmSq32.GetValue(numu_calc);
544 
545  double pre_seed_th23_LO = ( pre_seed_th23 > maxmixing ) ?
546  (2*maxmixing - pre_seed_th23) : pre_seed_th23;
547  double pre_seed_th23_UO = ( pre_seed_th23 > maxmixing ) ?
548  pre_seed_th23 : (2*maxmixing - pre_seed_th23);
549  // some fits were landing between 0.5 and 0.514, putting our
550  // pre_seed_th23_LO seed in the upper octant (> .5). If this happens,
551  // fall back to 0.45
552  if(pre_seed_th23_LO > .5) pre_seed_th23_LO = .45;
553 
554  std::cout << "Pre seeds--> ssth23_LO = " << pre_seed_th23_LO
555  << "\t dmsq_NH = " << fabs(pre_seed_dmsq) << endl;
556  // done getting pre seeds
557  //-------------------------------------------------------------
558 
559  // Set up a fitter to find the best fits constrained to true value
560  MinuitFitter fittrue(&expt, trueVars, systs);
561 
562  // Initialize to arbitrary values
563  double tb_delta = 999;
564  double tb_ssth23 = 999;
565  double tb_th13 = 999;
566  double tb_dmsq = 999;
567  std::vector<double> tf_delta_seeds = {0., 0.5, 1., 1.5};
568  std::vector<double> tf_th23_seeds;
569  if(upperOctant) tf_th23_seeds = {pre_seed_th23_UO, 0.501};
570  else tf_th23_seeds = {pre_seed_th23_LO, 0.499};
571 
572  std::map<const IFitVar*, std::vector<double>> tf_seedPts;
573  if(!deltaSlice) tf_seedPts.insert(std::pair<const IFitVar*,
574  std::vector<double>>
575  (&kFitDeltaInPiUnits, tf_delta_seeds));
576  if(!ssth23Slice) tf_seedPts.insert(std::pair<const IFitVar*,
577  std::vector<double>>
578  (th23_var, tf_th23_seeds));
579  if(!th13Slice) tf_seedPts.insert(std::pair<const IFitVar*,
580  std::vector<double>>
581  ( &kFitSinSq2Theta13, {0.082}));
582  if(!dmsqSlice) tf_seedPts.insert(std::pair<const IFitVar*,
583  std::vector<double>>
584  (dmsq_var,
585  { nh? fabs(pre_seed_dmsq):
586  -1*fabs(pre_seed_dmsq)}));
587 
588  // Run the true fit
589  std::cout << "Starting True fit ----->" <<endl;
590  auxShifts.ResetToNominal();
591  double trueChi = fittrue.Fit(calc,
592  auxShifts,
593  tf_seedPts,
594  thisSystShift,
595  IFitter::kQuiet)->EvalMetricVal();
596 
597  // save best parameters for posterity
598  tb_delta = kFitDeltaInPiUnits.GetValue(calc);
599  tb_ssth23 = th23_var->GetValue(calc);
600  tb_th13 = kFitSinSq2Theta13.GetValue(calc);
601  tb_dmsq = dmsq_var->GetValue(calc);
602  std::cout << "True fit parameters: ssth23 = " << tb_ssth23;
603  std::cout << ", delta = " << tb_delta;
604  std::cout << ", ss2th13 = " << tb_th13;
605  std::cout << ", dmsq = " << tb_dmsq << endl;
606  // make sure true vars don't somehow get changed from the true fit
607  if(deltaSlice) tb_delta = trueX;
608  if(ssth23Slice) tb_ssth23 = trueX;
609  if(dmsqSlice) tb_dmsq = trueX;
610 
611  // Set up fitter for calculating the global best fit
612  std::vector<const IFitVar*> bestVars;
613  bestVars.push_back(&kFitDeltaInPiUnits);
614  bestVars.push_back(&kFitDmSq32);
615  bestVars.push_back(&kFitSinSq2Theta13);
616  bestVars.push_back(&kFitSinSqTheta23);
617 
618  MinuitFitter fit(&expt, bestVars, systs);
619  std::vector <double> bf_dmsq_seeds = {fabs(pre_seed_dmsq),
620  -1*fabs(pre_seed_dmsq)};
621  std::vector <double> bf_delta_seeds = {0., 0.5, 1.0, 1.5};
622  std::vector <double> bf_th23_seeds = {pre_seed_th23_UO,
623  maxmixing,
624  pre_seed_th23_LO};
625 
626  std::vector <double> bf_th13_seeds;
627  if(th13Slice)bf_th13_seeds.push_back(trueX);
628  else bf_th13_seeds.push_back(0.082);
629 
630  std::map <const IFitVar*, std::vector<double>> bf_seedPts =
631  {{&kFitDeltaInPiUnits, bf_delta_seeds},
632  {&kFitSinSqTheta23, bf_th23_seeds},
633  {&kFitSinSq2Theta13, bf_th13_seeds},
634  {&kFitDmSq32, bf_dmsq_seeds}};
635 
636 
637  // Run Best fitter
638  // Initialize to arbitrary values.
639  double bestdelta = 999;
640  double bestssth23 = 999;
641  double bestth13 = 999;
642  double bestdmsq = 999;
643  std::cout << "Starting best fit ---->" << endl;
644  ResetOscCalcToDefault(calc);
645  auxShifts.ResetToNominal();
646  double bestChi = fit.Fit(calc,
647  auxShifts,
648  bf_seedPts,
649  thisSystShift,
650  IFitter::kQuiet)->EvalMetricVal();
651 
652  bestdelta = kFitDeltaInPiUnits.GetValue(calc);
653  bestssth23 = kFitSinSqTheta23.GetValue(calc);
654  bestth13 = kFitSinSq2Theta13.GetValue(calc);
655  bestdmsq = kFitDmSq32.GetValue(calc);
656  std::cout << "Best fit parameters: ssth23 = " << bestssth23;
657  std::cout << ", delta = " << bestdelta;
658  std::cout << ", ss2th13 = " << bestth13;
659  std::cout << ", dmsq = " << bestdmsq << endl;
660 
661 
662  // Make pseudoexperiment dst
663  FCPoint pt(seeddelta, seedssth23, seedss2th13, seeddmsq32,
664  pre_seed_th23, pre_seed_dmsq,
665  tb_delta, tb_ssth23, tb_th13, tb_dmsq,
666  trueChi,
667  bestdelta, bestssth23, bestth13, bestdmsq,
668  bestChi);
669 
670  fccol.AddPoint(pt);
671  unsigned int nfc_pts = fccol.NPoints();
672 
673  if(npts_flush != 0) {
674 
675  if(nfc_pts % npts_flush == 0 ) {
676  std::cout << "Saving " + std::to_string(nfc_pts) + "th point to " << fcout << endl;
677  fccol.SaveToFile(fcout);
678  }
679  }
680  }
681 
682  // Save to file
683  // Do our very best to insure we never have a same-named file and over-write
684  // good FC points.
685  }// end loop over bins
686  fccol.SaveToFile(fcout);
687 
688 }
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void ResetOscCalcToDefault(osc::IOscCalcAdjustable *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:23
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
double GetValue(const osc::IOscCalcAdjustable *osc) const override
double th23
Definition: runWimpSim.h:98
double delta
Definition: runWimpSim.h:98
double GetValue(const osc::IOscCalcAdjustable *osc) const override
Forward to wrapped Var&#39;s GetValue()
const FitDmSq32 kFitDmSq32
Definition: FitVars.cxx:18
void make_fc_slices_nersc_2018(int NPts, int startidx, int endidx, bool nh, int N, int jid, int npts_flush=0, TString plot="delta_LO", TString beam="both", bool corrSysts=false, bool test=false)
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Definition: FitVars.cxx:177
double GetValue(const osc::IOscCalcAdjustable *osc) const override
Definition: FitVars.cxx:171
const FitSinSqTheta23UpperOctant kFitSinSqTheta23UpperOctant
Definition: FitVars.cxx:16
void AddPoint(const FCPoint &p)
Definition: FCCollection.h:17
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Forward to wrapped Var&#39;s SetValue()
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)
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:333
const NeutronVisEScaleSyst2018 kNeutronVisEScalePrimariesSyst2018(false)
static SystShifts Nominal()
Definition: SystShifts.h:34
double GetValue(const osc::IOscCalcAdjustable *osc) const override
Definition: FitVars.cxx:42
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
const double kAna2018RHCPOT
Definition: Exposures.h:208
unsigned int NPoints() const
Definition: FCCollection.h:26
Represents the results of a single Feldman-Cousins pseudo-experiment.
Definition: FCPoint.h:8
const XML_Char * s
Definition: expat.h:262
std::vector< const ISyst * > getAllAna2018Systs(const EAnaType2018 ana, const bool smallgenie, const BeamType2018 beam, bool isFit)
const DummyAnaSyst kAnaCherenkovSyst("Cherenkov","Cherenkov")
expt
Definition: demo5.py:34
void SetSeeds(osc::IOscCalcAdjustable *calc, double delta, double th23, double dmsq)
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
std::string getenv(std::string const &name)
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:64
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
const double j
Definition: BetheBloch.cxx:29
float bin[41]
Definition: plottest35.C:14
std::vector< std::pair< TH1D *, double > > GetNumuCosmics2018(const int nq=4, std::string beam="fhc", bool GetFromUPS=false, bool NERSC=false)
OStream cout
Definition: OStream.cxx:6
const Binning bins
Definition: NumuCC_CPiBin.h:8
Combine multiple component experiments.
::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:384
std::pair< TH1D *, double > GetNueCosmics2018(std::string beam, bool GetFromUPS=false, bool NERSC=false)
virtual void SetValue(osc::IOscCalcAdjustable *osc, double val) const =0
const FitSinSqTheta23LowerOctant kFitSinSqTheta23LowerOctant
Definition: FitVars.cxx:17
void ResetToNominal()
Definition: SystShifts.cxx:141
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
assert(nhit_max >=nhit_nbins)
Interface definition for fittable variables.
Definition: IFitVar.h:16
std::vector< std::pair< const ISyst *, const ISyst * > > GetCorrelations(bool isNue, bool isFHC, bool ptExtrap)
const FitSinSq2Theta13 kFitSinSq2Theta13
Definition: FitVars.cxx:13
const IPrediction * GetNuePrediction2018(std::string decomp, osc::IOscCalc *calc, bool corrSysts, std::string beam, bool isFakeData, bool GetFromUPS=false, bool minimizeMemory=false, bool NERSC=false)
const double kAna2018FHCPOT
Definition: Exposures.h:207
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
void SaveToFile(const std::string &fname) const
const FitDmSq32InvertedHierarchy kFitDmSq32InvertedHierarchy
Definition: FitVars.cxx:23
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
const FitDmSq32NormalHierarchy kFitDmSq32NormalHierarchy
Definition: FitVars.cxx:20
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:78
const double kAna2018FHCLivetime
Definition: Exposures.h:213
void SetSystCorrelations(int idx, const std::vector< std::pair< const ISyst *, const ISyst * >> &corrs)
const double kAna2018RHCLivetime
Definition: Exposures.h:214
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)
Collection of FCPoint. Serializable to/from a file.
Definition: FCCollection.h:12
virtual double GetValue(const osc::IOscCalcAdjustable *osc) const =0
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17