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