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