make_fc_surfaces_nersc_2018.C
Go to the documentation of this file.
1 //#include <iostream>
2 #include <fstream>
3 
5 #include "CAFAna/Core/Spectrum.h"
7 
10 
13 
14 #include "CAFAna/Vars/FitVars.h"
15 
16 #include "CAFAna/FC/FCPoint.h"
17 #include "CAFAna/FC/FCCollection.h"
18 #include "CAFAna/FC/FCSurface.h"
19 
23 
25 
26 #include "CAFAna/Analysis/Calcs.h"
27 #include "OscLib/IOscCalc.h"
28 
29 #include "TRandom3.h"
30 #include "TH1.h"
31 
32 using namespace ana;
33 // Testing commit
34 // Needed for running on the OSG in a tag prior to r20638
35 
36 void SetSeeds(osc::IOscCalcAdjustable* calc, double delta, double th23, double dmsq)
37 {
39  kFitDeltaInPiUnits.SetValue(calc, delta);
40  kFitSinSqTheta23.SetValue(calc, th23);
41  kFitDmSq32.SetValue(calc, dmsq);
42 }
43 
44 void make_fc_surfaces_nersc_2018(int NPts, int startidx, int endidx, bool nh,
45  int N, int jid, int npts_flush = 0,
46  std::string plot = "deltassth23" ,
47  TString beam = "both", bool corrSysts = false,
48  bool readFromFile = false)
49 {
50  assert(plot == "deltath13" || plot=="ssth23dmsq32"
51  || plot == "deltassth23" && "Plot must be ssth23dmsq32 or deltassth23");
52 
53  assert(beam == "FHCOnly" || beam == "RHCOnly" || beam == "both");
54 
55  if(readFromFile && !corrSysts){
56  std::cout << "No bin list for stats only contours" << endl;
57  // exit(1);
58  }
59 
60  bool RHCOnly = beam.Contains("RHCOnly");
61  bool FHCOnly = beam.Contains("FHCOnly");
62  bool both = beam.Contains("both");
63  // ssth23dmsq32 has 400 total bins while deltassth23 has 900
64  int TotBins = 0;
65  TotBins = (plot == "deltassth23") ? 900 : 400;
66  int bins[TotBins];
67  std::memset(bins, 0, TotBins);
68  std::string hier = nh? "_NH_" : "_IH_";
69  std::string stat = corrSysts? "_systs": "_stats";
70 
71 
72  // If we don't want to run full corrections, we'll read from a text file which holds
73  // a list of bins to be corrected based on deltachisq in each bin. This list is
74  // generated by generate_fc_binlists.C
75 
76  ifstream binlist;
77  std::string path = getenv("FCHELPERANA2018_LIB_PATH");
78  path += "/binlists/";
79  TString binlist_fname = path+plot+hier+beam+stat+"_binlist.txt";
80  binlist.open(binlist_fname);
81  int j = 0;
82  // check if the binlist exists
83  if(!binlist.good()){
84  std::cout << "Could not find "<< binlist_fname << endl;
85  exit(1);
86  }
87  if(readFromFile){
88  do
89  {
90  binlist >> bins[j];
91  j++;
92  }while(!binlist.eof());
93  binlist.close();
94  }
95  else for (int i = 0; i < TotBins; i++) bins[i]=i;
96 
97  int ibin = bins[startidx];
98  int fbin = bins[endidx];
99 
100  std::string binname = "bins_"+std::to_string(ibin) + "_to_" + std::to_string(fbin);
101  std::string number = std::to_string(N);
102  std::string fcout = "FCCol_"+plot+"_"+hier+"_";
103  if(both) fcout+="both";
104  else fcout += RHCOnly? "RHCOnly":"FHCOnly";
105  fcout += corrSysts? "_systs_":"_stats_";
106  fcout += binname+"_v"+number+'_'+std::to_string(jid)+".root";
107 
108  std::string hierStr = nh ? "NH" : "IH";
109 
110  FCCollection fccol;
112 
113  // Relevant exposures, to be used later
114  double thisPOT;
115  double thisLivetime;
116 
117  // bools to make the function calls a little more readable
118  bool nueOnly = true;
119  bool numuOnly = true;
120  bool isFHC = true;
121  bool isRHC = true;
122  bool NERSC = true;
123 
124  // Grab the systs
125  std::vector<const ISyst*> systs = {};
126  if(corrSysts) systs = GetJointFitSystematicList(true, !nueOnly,
127  !numuOnly, isFHC, isRHC);
128  // need to grab different systs for nue rhc and fhc
129  std::vector<const ISyst*> snue_fhc = {};
130  if(corrSysts) snue_fhc = getAllAna2018Systs(kNueAna2018, true, kFHC);
131  std::vector<const ISyst*> snue_rhc = {};
132  if(corrSysts) snue_rhc = getAllAna2018Systs(kNueAna2018, true, kRHC);
133 
134  // numu systs are independent of beam mode
135  std::vector<const ISyst*> snumu = {};
136  if(corrSysts) snumu = GetJointFitSystematicList(true, !nueOnly,
137  numuOnly,
138  isFHC,
139  isRHC);
140 
141  std::vector <const IPrediction*> nue_preds;
142  std::vector <std::pair <TH1D*, double > > nue_cosmics;
143  std::vector <Spectrum> nue_cosmics_spectra;
144 
145 
146  // We're setting up the requisite inputs to the nue experiment
147 
148  //////////////////////////////////////////////////////////////
149  // Nue Predictions
150  /////////////////////////////////////////////////////////////
151  if(!RHCOnly)
152  nue_preds.push_back(GetNuePrediction2018("combo", calc, corrSysts, "fhc",
153  false, false, true, true));
154  if(!FHCOnly)
155  nue_preds.push_back(GetNuePrediction2018("prop", calc, corrSysts, "rhc",
156  false, false, true, true));
157  std::cout << "\n------------------Done Getting nue preds ---------------------\n";
158  if(!RHCOnly){
159  nue_cosmics.push_back(GetNueCosmics2018("fhc", false, true));
160  nue_cosmics_spectra.push_back(Spectrum(nue_cosmics.back().first,
162  }
163  if(!FHCOnly) {
164  nue_cosmics.push_back(GetNueCosmics2018("rhc", false, true));
165  nue_cosmics_spectra.push_back(Spectrum(nue_cosmics.back().first,
167  }
168 
169  ///////////////////////////////////////////////////////////////
170  // Numu predictions
171  //////////////////////////////////////////////////////////////
172 
173  // Construct numu inputs
174  std::vector < const IPrediction*> numu_preds;
175  std::vector < std::pair < TH1D*, double > > numu_cosmics;
176  std::vector <Spectrum> numu_cosmics_spectra;
177 
178  std::vector <const IPrediction*> this_numu_preds;
179  std::vector <std::pair<TH1D*, double > > this_numu_cosmics;
180 
181  std::vector<std::string> beams;
182  if(!RHCOnly) beams.push_back("fhc");
183  if(!FHCOnly) beams.push_back("rhc");
184  // Collect the eight numu preds/cosmics, 4 per beam mode
185  for(std::string beam : beams ) {
186  this_numu_preds = GetNumuPredictions2018(4, corrSysts, beam, false, kNuMu, true, true);
187  this_numu_cosmics = GetNumuCosmics2018(4, beam, false, true);
188  numu_preds.insert(numu_preds.end(),
189  this_numu_preds.begin(),
190  this_numu_preds.end());
191  numu_cosmics.insert(numu_cosmics.end(),
192  this_numu_cosmics.begin(),
193  this_numu_cosmics.end());
194  // Again, FHC and RHC have different pot and livetimes so
195  // cosmic hists need to be scaled appropriately
196  for(std::pair <TH1D*, double> h : this_numu_cosmics) {
197  thisPOT = kAna2018FHCPOT;
198  thisLivetime = kAna2018FHCLivetime;
199  if(beam == "rhc") {
200  thisPOT = kAna2018RHCPOT;
201  thisLivetime = kAna2018RHCLivetime;
202  }
203  numu_cosmics_spectra.push_back(Spectrum(h.first, thisPOT, thisLivetime));
204  }
205  }
206 
207  // Grab analysis contours to find bin center
208  std::string p = ""; //delta && !corrSysts
209  // if(plot == "ssth23dmsq32" && corrSysts) p = "ssth23dmsq32_";
210  // if(plot == "deltassth23" && corrSysts) p = "deltassth23_";
211  if(plot == "ssth23dmsq32") p = "dmsq_";
212  if(plot == "deltath13") p = "th13_";
213 
214  std::string helperupsname = std::string(getenv("FCHELPERANA2018_LIB_PATH"));
215  TString anasurfname = helperupsname
216  + "/contours/hist_contours_2018_joint_realData_"+beam+"_only";
217  anasurfname += nh? "NHcombo_" : "IHcombo_";
218  anasurfname += p;
219  anasurfname += corrSysts? "systs.root": "stats.root";
220 
221  TFile *anasurffile = TFile::Open(anasurfname);
222  anasurffile->cd();
223  std::string plotShortName="delta_th23";
224  if(plot == "ssth23dmsq32") plotShortName = "th23_dm32";
225  if(nh) plotShortName += "_NH";
226  else plotShortName += "_IH";
227  TH2* anasurf = (TH2F*)anasurffile->Get(plotShortName.c_str());
228 
229  // Grab the FCHelper2018 file
230  //std::string helperupsname = std::string(getenv("FCHELPERANA2017_LIB_PATH"));
231  TString helpername =
232  helperupsname+"FCInput/contours_FCInput_2018_joint_realData_"+beam+"_only";
233  helpername += nh? "NHcombo_" : "IHcombo_";
234  helpername += p;
235  helpername += corrSysts? "systs.root":"stats.root";
236 
237  TFile *fchelp = TFile::Open(helpername);
238 
239  // Find the bin center for this job, probably should be arguments
240  // dCP / pi
241 
242  // The surfaces are made with ana::ExpandedHistrogram() which sets the bin edges so that
243  // the input min/max x and y fall in the center of the bins. Grab these boundaries.
244  int NX = anasurf->GetXaxis()->GetNbins();
245  int NY = anasurf->GetYaxis()->GetNbins();
246  double minX = anasurf->GetXaxis()->GetBinCenter(1); // root bins start at 1
247  double maxX = anasurf->GetXaxis()->GetBinCenter(NX);
248  double minY = anasurf->GetYaxis()->GetBinCenter(1);
249  double maxY = anasurf->GetYaxis()->GetBinCenter(NY);
250 
251  double widthX = anasurf->GetXaxis()->GetBinWidth(1);
252  double widthY = anasurf->GetYaxis()->GetBinWidth(1);
253  std::cout << "widthX = " << widthX << endl;
254  std::cout << "widthY = " << widthY << endl;
255  assert(fbin < NX*NY && "Invalid bin number");
256 
257  int bin = 0;
258  for (int i = startidx; i <= endidx; i++) {
259  bin = bins[i];
260  std::cout<<"- "<<i<<std::endl;
261  std::cout << "bin " << bin << endl;
262  // Once we run out of important bins, get out.
263  if(bins[i] == 0 && i != 0) continue;
264 
265  // bin = bins[i];
266 
267  int binY = bin/NX;
268  int binX = bin - NX*binY;
269  std::cout << "binX: " << binX << endl;
270  std::cout << "binY: " << binY << endl;
271  anasurffile->cd();
272  double centerX = anasurf->GetXaxis()->GetBinCenter(binX+1);
273  double centerY = anasurf->GetYaxis()->GetBinCenter(binY+1);
274  std::cout << "centerX = " << centerX << endl;
275  std::cout << "centerY = " << centerY << endl;
276  fchelp->cd();
277 
278  // Pull out hists and get values at current bin
279  double seedssth23 = 0;
280  if(plot == "deltath13"){
281  TH2* hssth23 = (TH2F*) fchelp->Get((hierStr+"_SinSqTheta23").c_str());
282  seedssth23 = hssth23->GetBinContent(binX+1,binY+1);
283  delete hssth23;
284  }
285  double seeddmsq32 = 0;
286  if (plot == "deltassth23"){
287  TH2* hdmsq = (TH2F*)fchelp->Get((hierStr+"_DmSq32").c_str());
288  seeddmsq32 = hdmsq->GetBinContent(binX+1,binY+1); // root bins off-by-one
289  delete hdmsq;
290  }
291  double seeddelta = 0;
292  if (plot == "ssth23dmsq32"){
293  TH2* hdelta = (TH2F*)fchelp->Get((hierStr+"_DeltaCPpi").c_str());
294  seeddelta = hdelta->GetBinContent(binX+1,binY+1);
295  delete hdelta;
296  }
297  TH2* hss2th13 = (TH2F*)fchelp->Get((hierStr+"_SinSq2Theta13").c_str());
298  double seedss2th13 = hss2th13->GetBinContent(binX+1,binY+1);
299  delete hss2th13;
300 
301  // Get current syst shifts
302  std::map<std::string,double> seedsyst;
303  for (const ISyst* syst : systs){
304  TH2* h = (TH2F*)fchelp->Get((hierStr+"_"+syst->ShortName()).c_str());
305  if (!h){
306  std::cout << "Don't see a prof hist for " << syst->ShortName() << ". Continuing, but check your ups version." << std::endl;
307  continue;
308  }
309  double seedval = h->GetBinContent(binX+1,binY+1);
310  seedsyst.insert(std::pair<std::string,double>(syst->ShortName(),seedval));
311  }
312 
313  std::cout << "Parameters from anasurf " << endl;
314  std::cout << "ssth23 = " << seedssth23 << endl;
315  std::cout << "dmsq = " << seeddmsq32 << endl;
316  std::cout << "delta = " << seeddelta << endl;
317  std::cout << "th13 = " << seedss2th13 << endl;
318 
319  // Need random numbers to sample the bin
320  TRandom3 rnd(0);
321 
322  for (int i = 0; i < NPts; i++){
323  ResetOscCalcToDefault(calc);
324 
325  const IFitVar* dmsq_var;
326  if(nh) dmsq_var = &kFitDmSq32ScaledNH;
327  else dmsq_var = &kFitDmSq32ScaledIH;
328 
329  // backwards compatability
330  double trueX = centerX;
331  double trueY = centerY;
332  std::cout << "trueX: " << trueX << endl;
333  std::cout << "trueY: " << trueY << endl;
334  // Plot = deltassth23: x=delta, y=ssth23
335  // Plot = ssth23dmsq: x=ssth23, y=dmsq32
336  // Plot = deltassth13: x=th13, y=delta
337  if(plot != "deltath13"){
338  kFitDeltaInPiUnits.SetValue(calc, plot=="deltassth23"?trueX:seeddelta);
339  kFitSinSqTheta23.SetValue(calc, plot=="deltassth23"?trueY:trueX);
340  dmsq_var->SetValue(calc, plot=="deltassth23"?seeddmsq32:trueY);
341  kFitSinSq2Theta13.SetValue(calc, seedss2th13);
342  }
343  else{
344  kFitDeltaInPiUnits.SetValue(calc, trueY);
345  kFitSinSqTheta23.SetValue(calc, seedssth23);
346  dmsq_var->SetValue(calc, seeddmsq32);
347  kFitSinSq2Theta13.SetValue(calc, trueX);
348  }
349 
350  // Need to seed psuedo-experiments at profiled syst shifts read from
351  // FCInput
352  SystShifts seednue_fhc;
353  for (const ISyst *syst : snue_fhc){
354  if (seedsyst.find(syst->ShortName()) == seedsyst.end()) continue;
355  seednue_fhc.SetShift(syst,seedsyst[syst->ShortName()]);
356  }
357  SystShifts seednue_rhc;
358  for (const ISyst *syst : snue_rhc) {
359  if (seedsyst.find(syst->ShortName()) == seedsyst.end()) continue;
360  seednue_rhc.SetShift(syst, seedsyst[syst->ShortName()]);
361 
362  }
363  SystShifts seednumu;
364  for (const ISyst *syst : snumu){
365  if (seedsyst.find(syst->ShortName()) == seedsyst.end()) continue;
366  seednumu.SetShift(syst,seedsyst[syst->ShortName()]);
367  }
368 
369 
370  // Make a mock data spectrum for nue
371  // Use mock data + cosmic as predictions
372 
373  std::cout << "----------Making mock data at-----------------\n";
374  std::cout << "ssth23 = " << kFitSinSqTheta23.GetValue(calc);
375  std::cout << ", delta = " << kFitDeltaInPiUnits.GetValue(calc);
376  std::cout << ", dmsq = " << dmsq_var->GetValue(calc);
377  std::cout << ", th13 = " << kFitSinSq2Theta13.GetValue(calc);
378  std::cout << endl;
379 
380  std::vector<Spectrum> mocknue;
381  SystShifts this_seednue;
382  int mock_integral = 0;
383  for(int j = 0; j < (int) nue_preds.size(); j++) {
384  if(j == 0 && !RHCOnly){
385  this_seednue = seednue_fhc;
386  thisPOT = kAna2018FHCPOT;
387  }
388  else{
389  this_seednue = seednue_rhc;
390  thisPOT = kAna2018RHCPOT;
391  }
392  Spectrum fakenue = nue_preds[j]->Predict(calc);
393  if(corrSysts) fakenue = nue_preds[j]->PredictSyst(calc, this_seednue);
394  fakenue += nue_cosmics_spectra[j];
395  mocknue.push_back(fakenue.MockData(thisPOT));
396  // Save this for posterity's sake
397  mock_integral += mocknue[j].Integral(thisPOT);
398  }
399  std::cout << "------------Made "
400  << nue_preds.size()
401  << " mock nue data spectra ---------------\n";
402 
403  // Make a mock data spectrum for numu
404  std::vector<Spectrum> mocknumu;
405  for (int j = 0; j < (int) numu_preds.size(); j++){
406 
407  if(j < 3 && !RHCOnly)
408  thisPOT = kAna2018FHCPOT;
409  else
410  thisPOT = kAna2018RHCPOT;
411 
412  Spectrum fakenm = numu_preds[j]->Predict(calc);
413  if(corrSysts) numu_preds[j]->PredictSyst(calc,seednumu);
414  fakenm += numu_cosmics_spectra[j];
415  mocknumu.push_back(fakenm.MockData(thisPOT));
416  }
417 
418  std::cout << "------------Made "
419  << numu_preds.size()
420  << " mock numu data spectra ----------------\n";
421 
422  // Set up experiment, nova nue + nova numu + 2017 reactor
423  // Get true chi-square before fit
424 
425  //TODO experiments
426  std::vector<const IExperiment*> expts;
427  std::vector<const IExperiment*> numu_expts;
428  for(int j = 0; j < (int)nue_preds.size(); j++) {
429  expts.push_back(new SingleSampleExperiment(nue_preds[j], mocknue[j],
430  nue_cosmics[j].first,
431  nue_cosmics[j].second)); // Nue
432  }
433  for (int j = 0; j < (int)numu_preds.size(); j++) {
434  expts.push_back(new SingleSampleExperiment(numu_preds[j], mocknumu[j],
435  numu_cosmics[j].first,
436  numu_cosmics[j].second));
437  numu_expts.push_back(new SingleSampleExperiment(numu_preds[j], mocknumu[j],
438  numu_cosmics[j].first,
439  numu_cosmics[j].second));
440  }
441  if(plot != "deltath13") {
442  expts.push_back( WorldReactorConstraint2017() );
443  numu_expts.push_back( WorldReactorConstraint2017() );
444  }
445  std::cout << "Creating multiexperiment from " << expts.size() << " experiments" << endl;
446  MultiExperiment expt(expts);
447  MultiExperiment numu_expt(numu_expts);
448  // Correlate your systs
449  // TODO GetCorrelation function calls are different
450  for(int j = 0 ; j < (int)nue_preds.size(); j++){
451  if(j==0 && !RHCOnly)
452  expt.SetSystCorrelations(j, GetCorrelations(true, true)); // nue FHC
453  else
454  expt.SetSystCorrelations(j, GetCorrelations(true, false)); // nue RHC
455  auto notfornumu = GetCorrelations(false, false); // numu not dependent on beam
456  for(int j = 0; j < (int)numu_preds.size(); ++j)
457  expt.SetSystCorrelations(j+nue_preds.size(), notfornumu);
458  }
459 
460 
461  /////////////////////////////////////////////////////////////////
462  // Numu only fit to figure out dmsq seeds
463  // (may be able to do without this for contours)
464  /////////////////////////////////////////////////////////////////
465  std::cout << "Fitting numu to find pre seeds" << endl;
466 
467  double maxmixing = .514;
469  SystShifts auxShifts = SystShifts::Nominal();
470 
471  MinuitFitter fitnumu_only(&numu_expt, {&kFitSinSqTheta23, &kFitDmSq32});
472 
473  double minchi_numu = fitnumu_only.Fit(numu_calc, auxShifts,
474  {{&kFitSinSqTheta23, {0.4}},
475  {&kFitDmSq32, {nh? 2.35e-3 : -2.35e-3}}},
476  {},
477  IFitter::kQuiet)->EvalMetricVal();
478 
479  double pre_seed_dmsq = dmsq_var->GetValue(numu_calc);
480  double pre_seed_th23 = kFitSinSqTheta23.GetValue(numu_calc);
481 
482  double pre_seed_th23_LO = ( pre_seed_th23 > maxmixing ) ?
483  (2*maxmixing - pre_seed_th23) : pre_seed_th23;
484  double pre_seed_th23_UO = ( pre_seed_th23 > maxmixing ) ?
485  pre_seed_th23 : (2*maxmixing - pre_seed_th23);
486  if(pre_seed_th23_LO > .5) pre_seed_th23_LO = .45;
487 
488  std::cout << "Preseeds ---> ssth23_LO = " << pre_seed_th23_LO;
489  std::cout << "\tdmsq_NH = " << fabs(pre_seed_dmsq) << endl;
490  // done getting pre seeds
491  // --------------------------------------------------------------
492  ////////////////////////////////////////////////////////////////////
493  // True fit
494  ////////////////////////////////////////////////////////////////////
495 
496  std::vector<const IFitVar*> trueVars;
497  if (plot=="deltassth23"){
498  trueVars.push_back(dmsq_var);
499  trueVars.push_back(&kFitSinSq2Theta13);
500  }
501  if (plot=="ssth23dmsq32"){
502  trueVars.push_back(&kFitDeltaInPiUnits);
503  trueVars.push_back(&kFitSinSq2Theta13);
504  }
505 
506 
507  MinuitFitter fittrue(&expt, trueVars, systs);
508  double tb_th13 = 0;
509  double tb_dmsq = 0;
510  double tb_ssth23 = 0;
511  double tb_delta = 0;
512  SystShifts systshift = {};
513  // if not correcting deltassth23, then put one element in this
514  // vector so the following loop will run
515  std::vector <double> tf_delta_seeds = {0, .5, 1., 1.5};
516  std::vector <double> tf_ssth23_seeds = {pre_seed_th23_LO,
517  maxmixing,
518  pre_seed_th23_UO};
519 
520  // make a map of all of the seeds for the true fit
521  std::map<const IFitVar*, std::vector<double>> tf_seedPts;
522  if(plot == "deltassth23")
523  tf_seedPts = {{dmsq_var, {nh? fabs(pre_seed_dmsq) : -1*fabs(pre_seed_dmsq)}},
524  {&kFitSinSq2Theta13, {0.082}}};
525  else if(plot == "ssth23dmsq32")
526  tf_seedPts = {{&kFitDeltaInPiUnits, tf_delta_seeds},
527  {&kFitSinSq2Theta13, {0.082}}};
528  else if(plot == "deltath13")
529  tf_seedPts = {{dmsq_var, {nh? fabs(pre_seed_dmsq) : -1*fabs(pre_seed_dmsq)}},
530  {&kFitSinSqTheta23, tf_ssth23_seeds}};
531 
532  std::cout << "True fit parameters (pre fit): ssth23 = " << kFitSinSqTheta23.GetValue(calc);
533  std::cout << ", delta = " << kFitDeltaInPiUnits.GetValue(calc);
534  std::cout << ", ss2th13 = " << kFitSinSq2Theta13.GetValue(calc);
535  std::cout << ", dmsq = " << dmsq_var->GetValue(calc) << " x10^-3 eV" << endl;
536 
537 
538  std::cout << "Starting true fit --->" << endl;
539  auxShifts.ResetToNominal();
540  double trueChi = fittrue.Fit(calc,
541  auxShifts,
542  tf_seedPts,
543  {},
544  IFitter::kQuiet)->EvalMetricVal();
545 
546  // save best parameters for posterity
547  tb_delta = kFitDeltaInPiUnits.GetValue(calc);
548  tb_ssth23 = kFitSinSqTheta23.GetValue(calc);
549  tb_dmsq = dmsq_var->GetValue(calc);
550  tb_th13 = kFitSinSq2Theta13.GetValue(calc);
551 
552  std::cout << "True fit parameters (post fit): ssth23 = " << tb_ssth23;
553  std::cout << ", delta = " << tb_delta;
554  std::cout << ", ss2th13 = " << tb_th13;
555  std::cout << ", dmsq = " << tb_dmsq << " x10^-3 eV" << endl;
556 
557  // Set up best fitter
558  MinuitFitter fit(&expt,
561  systs);
562 
563  std::vector <double> bf_dmsq_seeds = {fabs(pre_seed_dmsq),
564  -1*fabs(pre_seed_dmsq)};
565  std::vector <double> bf_delta_seeds = {0., 0.5, 1., 1.5};
566  std::vector <double> bf_th23_seeds = {pre_seed_th23_LO, pre_seed_th23_UO,
567  maxmixing};
568  std::vector <double> bf_th13_seeds = {0.082};
569  std::map <const IFitVar*, std::vector<double>> bf_seedPts =
570  {{&kFitDeltaInPiUnits, bf_delta_seeds},
571  {&kFitSinSqTheta23, bf_th23_seeds},
572  {&kFitSinSq2Theta13, bf_th13_seeds},
573  {&kFitDmSq32Scaled, bf_dmsq_seeds}};
574 
575  double bestdelta = 0;
576  double bestssth23 = 0;
577  double bestth13 = 0;
578  double bestdmsq = 0;
579  std::cout << "Starting best fit ----> " << endl;
580  ResetOscCalcToDefault(calc);
581  auxShifts.ResetToNominal();
582  double bestChi = fit.Fit(calc,
583  auxShifts,
584  bf_seedPts,
585  {},
586  IFitter::kQuiet)->EvalMetricVal();
587 
588  bestdelta = kFitDeltaInPiUnits.GetValue(calc);
589  bestssth23 = kFitSinSqTheta23.GetValue(calc);
590  bestth13 = kFitSinSq2Theta13.GetValue(calc);
591  bestdmsq = kFitDmSq32Scaled.GetValue(calc);
592  std::cout << "Best fit parameters: ssth23 = " << bestssth23;
593  std::cout << ", delta = " << bestdelta;
594  std::cout << ", ss2th13 = " << bestth13;
595  std::cout << ", dmsq = " << bestdmsq << " x10^-3 eV"<< endl;
596 
597 
598  // Make pseudoexperiment dst
599  // with FCPoint in its current state (r31771) the names of the branches are
600  // not correct. We want to save best fit information
601  FCPoint pt(seeddelta, seedssth23, seedss2th13, seeddmsq32,
602  pre_seed_th23, pre_seed_dmsq,
603  tb_delta, tb_ssth23, tb_th13, tb_dmsq,
604  trueChi,
605  bestdelta, bestssth23, bestth13, bestdmsq,
606  bestChi);
607 
608  fccol.AddPoint(pt);
609 
610  unsigned int nfc_pts = fccol.NPoints();
611  if(npts_flush != 0) {
612  if(nfc_pts % npts_flush == 0) {
613  std::cout << "Saving " + std::to_string(nfc_pts) + "th point to " << fcout << endl;
614  fccol.SaveToFile(fcout);
615  }
616  }
617  }
618  }// end loop over bins
619  // Save to file
620  fccol.SaveToFile(fcout);
621 }
double minY
Definition: plot_hist.C:9
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
const FitDmSq32ScaledNH kFitDmSq32ScaledNH
Definition: FitVars.cxx:21
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
const char * p
Definition: xmltok.h:285
void AddPoint(const FCPoint &p)
Definition: FCCollection.h:17
double maxY
Definition: plot_hist.C:10
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:249
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Definition: FitVars.cxx:48
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
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
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
std::vector< const ISyst * > getAllAna2018Systs(const EAnaType2018 ana, const bool smallgenie, const BeamType2018 beam, bool isFit)
expt
Definition: demo5.py:34
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:69
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
const double j
Definition: BetheBloch.cxx:29
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
float bin[41]
Definition: plottest35.C:14
std::vector< float > Spectrum
Definition: Constants.h:610
std::vector< std::pair< TH1D *, double > > GetNumuCosmics2018(const int nq=4, std::string beam="fhc", bool GetFromUPS=false, bool NERSC=false)
static bool isFHC
OStream cout
Definition: OStream.cxx:6
const std::string path
Definition: plot_BEN.C:43
const Binning bins
Definition: NumuCC_CPiBin.h:8
Combine multiple component experiments.
const FitDmSq32ScaledIH kFitDmSq32ScaledIH
Definition: FitVars.cxx:22
std::pair< TH1D *, double > GetNueCosmics2018(std::string beam, bool GetFromUPS=false, bool NERSC=false)
virtual void SetValue(osc::IOscCalcAdjustable *osc, double val) const =0
void ResetToNominal()
Definition: SystShifts.cxx:143
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
exit(0)
assert(nhit_max >=nhit_nbins)
Interface definition for fittable variables.
Definition: IFitVar.h:16
void SetSeeds(osc::IOscCalcAdjustable *calc, double delta, double th23, double dmsq)
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
void make_fc_surfaces_nersc_2018(int NPts, int startidx, int endidx, bool nh, int N, int jid, int npts_flush=0, std::string plot="deltassth23", TString beam="both", bool corrSysts=false, bool readFromFile=false)
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
std::vector< const ISyst * > GetJointFitSystematicList(bool corrSysts, bool nueExclusive=false, bool numuExclusive=false, bool isFHC=true, bool isRHC=true, bool intersection=true, bool ptExtrap=true)
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:80
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
enum BeamMode string